function [eqdskval]=readeqdsk_simple(fname,COCOS_in);
%
% Read eqdsk files
%
% function [eqdskval]=readeqdsk(fname,varargin);
%
% if fname=file_name is empty or not given, prompts for it
%
% COCOS_in: COCOS index value (see paper in http://crpp.epfl.ch/COCOS for details). Assume 2 (as CHEASE) if not given
%

cocos_in = 2;
if nargin>1 && ~isempty(COCOS_in})
  cocos_in = COCOS_in;
end

% check file name and input
if nargin == 0
  [fname,pathname] = uigetfile([{'EQDSK*;eqdsk*'},{'EQDSK* or eqdsk* files'}; ...
                    {'*EQDSK*;*eqdsk*'},{'*EQDSK* or *eqdsk* files'};{'*'},{'All files'}],'Select an EQDSK file');
  if isequal(fname,0) || isequal(pathname,0); return; end
  fnamefull = fullfile(pathname, fname);
else;
  fnamefull = fname;
  if ~exist(fnamefull); 
    disp([fnamefull,' does not exist, prompts for it']); 
    [fname,pathname] = uigetfile([{'EQDSK*;eqdsk*'},{'EQDSK* or eqdsk* files'}; ...
                    {'*EQDSK*;*eqdsk*'},{'*EQDSK* or *eqdsk* files'};{'*'},{'All files'}],'Select an EQDSK file');
    if isequal(fname,0) || isequal(pathname,0); return; end
    fnamefull = fullfile(pathname, fname);
  else
    ii=findstr(fnamefull,filesep);
    if ~isempty(ii)
      fname=fnamefull(ii(end)+1:end);
      pathname=fnamefull(1:ii(end)-1);
    else
      fname=fnamefull;
      pathname=pwd;
    end
  end
end

eqdskval.fnamefull=fnamefull;
eqdskval.fname=fname;
eqdskval.pathname=pathname;
eqdskval.cocos = cocos_in;

% Get assumed sign relations given the COCOS index
cocos_struct = cocos(cocos_in);

fid=fopen(fnamefull,'r');
if fid==-1; error(['error opening ',fnamefull]); end

eqdskval.stitle=fscanf(fid,'%c',48);
ind=fscanf(fid,'%d',3);
% goes to next line aftar char(10)
nextchar=fscanf(fid,'%c',1);
ii=1;
while double(nextchar)~=10
  nextchar=fscanf(fid,'%c',1);
  ii=ii+1;
end
eqdskval.ind1=ind(1);
eqdskval.nr=ind(2);
eqdskval.nz=ind(3);
for i=1:4
  line(i,1:5)=fscanf(fid,'%g',5)';
end
eqdskval.F=fscanf(fid,'%g',eqdskval.nr);
eqdskval.p=fscanf(fid,'%g',eqdskval.nr);
eqdskval.FFprime=fscanf(fid,'%g',eqdskval.nr);
eqdskval.pprime=fscanf(fid,'%g',eqdskval.nr);
eqdskval.psirz=fscanf(fid,'%g',eqdskval.nr*eqdskval.nz);
eqdskval.q=fscanf(fid,'%g',eqdskval.nr);
eqdskval.psimesh=linspace(0,1,eqdskval.nr)';
eqdskval.rhopsi=sqrt(eqdskval.psimesh);
eqdskval.nbbound=fscanf(fid,'%d',1);
eqdskval.nblim=fscanf(fid,'%d',1);
rzplasma=fscanf(fid,'%g',2*eqdskval.nbbound);
eqdskval.rplas=rzplasma(1:2:end);
eqdskval.zplas=rzplasma(2:2:end);
rzlimiter=fscanf(fid,'%g',2*eqdskval.nblim);
eqdskval.rlim=rzlimiter(1:2:end);
eqdskval.zlim=rzlimiter(2:2:end);
%
% read extra lines as comments in extralines cell array
% footer of file
ii=1;
while 1
  tline = fgetl(fid);
  if ~ischar(tline);break;end;
  footer{ii}=tline;
  ii=ii+1;
end
if length(footer)>1 && isempty(footer{1})
  eqdskval.extralines = footer(2:end);
else
  eqdskval.extralines=footer;
end
% strvcat(footer)

fclose(fid);

%
eqdskval.psi=NaN*ones(eqdskval.nr,eqdskval.nz);
for j=1:eqdskval.nz
  eqdskval.psi(:,j)=eqdskval.psirz((j-1)*eqdskval.nr+[1:eqdskval.nr]);
end

eqdskval.rboxlen=line(1,1);
eqdskval.zboxlen=line(1,2);
eqdskval.r0=line(1,3);
eqdskval.rboxleft=line(1,4);
eqdskval.zmid=line(1,5);
eqdskval.raxis=line(2,1);
eqdskval.zaxis=line(2,2);
eqdskval.psiaxis=line(2,3);
eqdskval.psiedge=line(2,4);
eqdskval.b0=line(2,5);
eqdskval.ip=line(3,1);

eqdskval.rmesh=linspace(eqdskval.rboxleft,eqdskval.rboxleft+eqdskval.rboxlen,eqdskval.nr);
eqdskval.zmesh=linspace(eqdskval.zmid-eqdskval.zboxlen/2,eqdskval.zmid+eqdskval.zboxlen/2,eqdskval.nz);

% construct limiter from box if needed
if length(eqdskval.rlim)==1 & eqdskval.nblim==1
  eqdskval.rlim=[eqdskval.rboxleft ; eqdskval.rboxleft+eqdskval.rboxlen ; ...
    eqdskval.rboxleft+eqdskval.rboxlen ; eqdskval.rboxleft ; ...
    eqdskval.rboxleft];
  eqdskval.zlim=[-eqdskval.zboxlen/2 ; -eqdskval.zboxlen/2 ; ...
    +eqdskval.zboxlen/2 ; +eqdskval.zboxlen/2 ; ...
    -eqdskval.zboxlen/2];
  eqdskval.nblim = 5;
elseif length(eqdskval.rlim)==1
  disp(['Inconsistency: length(eqdskval.rlim)= ' num2str(length(eqdskval.rlim)) ...
        ' but eqdskval.nblim = ' num2str(eqdskval.nblim)]);
end

% check COCOS consistency
% check q sign
qsign=sign(eqdskval.q(end));
if qsign*cocos_struct.sigma_rhothetaphi*sign(eqdskval.ip)*sign(eqdskval.b0) < 0
  warning(['sign of q is: ' num2str(qsign)])
  warning(['But is should be equal to sign_rhothetaphi*sign(Ip)*sign(B0) = ' num2str(cocos_struct.sigma_rhothetaphi) '*' ...
        num2str(sign(eqdskval.ip)) '*' num2str(sign(eqdskval.b0)) ' = ' ...
        num2str(cocos_struct.sigma_rhothetaphi*sign(eqdskval.ip)*sign(eqdskval.b0))])
end
if sign(eqdskval.F)*sign(eqdskval.b0) < 0
  warning('Signs of F and B0 are not consistent')
end
if sign(eqdskval.psiedge-eqdskval.psiaxis)*cocos_struct.sigma_Bp*sign(eqdskval.ip) < 0
  if eqdskval.psiedge > eqdskval.psiaxis
    warning(['psi should be decreasing  with sign(Ip) = ' num2str(sign(eqdskval.ip)) ' for COCOS = ' num2str(cocos_struct.cocos)]);
  else
    warning(['psi should be increasing  with sign(Ip) = ' num2str(sign(eqdskval.ip)) ' for COCOS = ' num2str(cocos_struct.cocos)]);
  end
elseif sign(mean(eqdskval.pprime))*sign(eqdskval.ip)*cocos_struct.sigma_Bp > 0
  warning(['sign(pprime) should be ' num2str(-sign(eqdskval.ip)*cocos_struct.sigma_Bp)])
end

