% minimalist demo to run chease with matlab scripts
clear; close all;

%addpath /home/ffelici/matlab/CHEASEgui_develop

namelist_file = 'namelist_ITERhybrid';
EXPEQ_file = 'EXPEQ_ITERhybrid';

% load data from Bert Maljaars
load ITER_hybrid_citrin_equil
R0 = 6.2; B0 = 5.3;
RZcontours =  [eq_R(end,:)',eq_Z(end,:)'-eq_Z(1,1)]/R0;
psi = prof_psi-prof_psi(1); psi=psi/sign(psi(end));

psi_chease = psi/(R0^2*B0);
ptot_chease = ptot/(B0^2/(4e-7*pi));
F_chease = F/(R0*B0);

psiN = (psi-psi(1))/(psi(end)-psi(1));
rhopsi = sqrt(psiN);

rho_grid = linspace(0,1,201);
psi_grid = rho_grid.^2*psi_chease(end);
[T,Tprime] = interpos(psi_chease,F_chease,psi_grid,1e-12,[0 0],[0,0]);
[P,Pprime] = interpos(psi_chease,ptot_chease,psi_grid,1e-9,[1 2],[0 0]);
TTprime = T.*Tprime;

clf;
subplot(211);
plot(psi_chease,F_chease,'b',psi_grid,T,'r--')
subplot(212);
plot(psi_chease,ptot_chease,'b',psi_grid,P,'r--')

%

EXPEQ_data = write_EXPEQ; % default data
EXPEQ_data.RZ_psi = RZcontours;
EXPEQ_data.TTprime = TTprime;
EXPEQ_data.Pprime = Pprime;
EXPEQ_data.nsttp = 1;
EXPEQ_data.nrhotype = 0;
EXPEQ_data.rho = rho_grid;
%EXPEQ_data.Jparallel = [];
%EXPEQ_data.Jparallel = Jpar;

[EXPEQdataout,fclose_out] = write_EXPEQ(EXPEQ_data,EXPEQ_file);

[~,ha] = plot_EXPEQ(EXPEQdataout.fnamefull,'b.-'); % plot input EXPEQ file
%%
nl.nsttp = EXPEQ_data.nsttp; % set correct nsttp in namelist
nl.npropt = nl.nsttp; 
nl.b0exp = B0; nl.r0exp = R0;
nl.ncscal = 4; % no scaling, use input TT' and p'
nl = write_namelist_CHEASE(namelist_file,nl);

!PATH=$PATH:/home/sauter/bin
fname_out = run_chease(namelist_file,EXPEQ_file);

%%
plot_EXPEQ(EXPEQdataout.fnamefull,'b'); drawnow;

%
% find EXPEQ_OUT file
for ii=1:length(fname_out);
    if ~isempty(strfind(fname_out{ii},'EXPEQ.OUT_'))
        EXPEQ_OUT_file = fname_out{ii}; break
    end
end
% find cols file
for ii=1:length(fname_out);
    if ~isempty(strfind(fname_out{ii},'.cols'))
        ocolsfile = fname_out{ii}; break
    end
end
% find EQDSK
for ii=1:length(fname_out);
    if ~isempty(strfind(fname_out{ii},'EQDSK'))
        eqdskfile = fname_out{ii}; break
    end
end

% plot output EXPEQ file in same plot
plot_EXPEQ(EXPEQ_OUT_file,'r.-');
% it is normal that j|| is different

%
profiledata = read_ocols(ocolsfile);

% also read eqdsk data
eqdskdata = readeqdsk(eqdskfile);
cheasedata.eqdsk = eqdskdata;
cheasedata.profiledata = profiledata; 

save ITER_hybrid_citrin_equil_cheasedata.mat cheasedata

%%
% compare with CRONOS result
%plot(cheasedata.psichease_eq_psi_over_2pi.data*2*pi,cheasedata.qprofile.data,'b')

