%%% calculation of profiles after sawtooth crash - incomplete reconnection model 
%% % saw_incom module
function varargout = saw_incom(varargin)

if nargin == 0,
    module_params = struct(...
        's1crit', 0.2,... %critical shear
        'mode', 1, ...%reconnection model: 0 full, 1 incomplete
        'w_norm', 0.5,...
        'check', false...
        ); 
    varargout{1} = module_params;
    return %empty call, probably to get default structures
elseif nargin==7;
    x = varargin{1};
    %g = varargin{2};
    v = varargin{3};
    %u = varargin{4};
    %it = varargin{5};
    model  = varargin{6};
    saw_params = varargin{7}; % distribute inputs
else
    error('must call with 0 or 7 inputs');
end

xnew = x; 

Psihat = x(model.psi.xind);
Tehat = x(model.te.xind);

% %% workaround for reversed shear q
% % check negative shear
% iota = model.out.Ciotagauss*Psihat;
% iq1=find(iota>=1,1,'last'); % outermost point with iota >=1
% iqrev = find(iota(1:iq1)<1,1,'last'); % index inside this region where iota<1 again
% 
% if ~isempty(iqrev) % this means there is a reversed shear region
%     iotafix=iota;
%     iotafix(1:iqrev) = iota(iqrev+1);
%     Psihatnew = [model.out.Ciotagauss;model.psi.Lam(end,:)]\[iotafix;model.psi.Lam(end,:)*Psihat];
%     Psihat = Psihatnew;
% end

%% workaround for reversed shear q (similar to Kim's option 1 in ASTRA)
% % check negative shear
iota = model.out.Ciotagauss*Psihat;
iota_th = 1.002; % fixed iota value for iota < 1.
iq1 = find(iota >= 1., 1, 'last'); % outermost point with iota >=1
rev = (iota(1:iq1) < 1.); % iota < 1 gaps

if ~isempty(rev) % this means there is a reversed shear region
    iotafix = iota;
    iotafix(rev) = iota_th;
    Psihatnew = [model.out.Ciotagauss;model.psi.Lam(end,:)]\[iotafix;model.psi.Lam(end,:)*Psihat];
    Psihat = Psihatnew;
end

if saw_params.check
    hs=subplot(221);
    plot(hs,model.rgrid.rhogauss,1./iota,'b.-',...
        model.rgrid.rhogauss(iq1),1./iota(iq1),'rx',...
        model.rgrid.rhogauss(rev),1./iota(rev),'go'...
        );
    drawnow
end
%%


% metric
rhogauss = model.rgrid.rhogauss;

% finer grid
r = linspace(rhogauss(1), rhogauss(end), 1001);

% minor radius
rtor_gauss = model.profiles.rtor(end) / 60. *(1:60)';
% aminor_gauss = model.equi.epsilon * model.equi.R0 .* rhogauss; % minor radius
aminor_gauss = rhogauss .* rtor_gauss(end);
aminor = interp1(rhogauss,aminor_gauss,r,'cubic');

B0 = model.equi.B0;
psi_gauss = model.psi.Lamgauss * Psihat; 
% phi_gauss = pi .* B0 * (rhogauss .* rtor_gauss(end)) .^2;
phi_gauss = pi .* B0 * aminor_gauss .^2;

% helical flux
psis_gauss = (psi_gauss - psi_gauss(1)) - (phi_gauss - phi_gauss(1));     
psis = interp1(rhogauss, psis_gauss, r, 'cubic');

% volume
vol_gauss = model.profiles.Vgauss;
vol = interp1(rhogauss, vol_gauss, r, 'cubic');

%electron density
ne_gauss = model.te.LamTegauss * v(model.ne.vind);
ne = interp1(rhogauss, ne_gauss, r, 'cubic');

%electron temperature
Te_gauss = model.te.LamTegauss * Tehat;
Te = interp1(rhogauss, Te_gauss, r, 'cubic');


%% define regions: 
% 1: 0 to psistar max
% 2: psistar max to psistar=0 (mixing radius)

% define 1
[~,imax] = max(psis);
ii1 = 1:imax;
 
% define 2
imix = find((psis(1:end-1)>0) & (diff(psis)<0), 1, 'last') + 1;
ii2 = imax:imix;

%evaluate quantities at imix
r_mix = r(imix-1) + (0-psis(imix-1))*(r(imix)-r(imix-1))/(psis(imix)-psis(imix-1));
a_mix = interp1(r,aminor,r_mix);
vol_mix = interp1(r,vol,r_mix);
ne_mix = interp1(r, ne, r_mix);
Te_mix = interp1(r,Te,r_mix);
psis_mix = 0; % by definition

%add quantities at imix
p = imix-1;
r_i = insert_1point(r,r_mix,p);
a_i = insert_1point(aminor,a_mix,p);
vol_i = insert_1point(vol,vol_mix,p);
psis_i = insert_1point(psis,psis_mix,p);
ne_i = insert_1point(ne,ne_mix,p);
Te_i = insert_1point(Te,Te_mix,p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
w_crit = saw_params.w_norm * a_i(imix);% critical island width : r2-r1 = w_crit
a1 = a_i(ii1);

% find a2 corresponding to each psis(a1)
ppp = psis_i(imix+imax-(ii2)); %from imix to imax in reversed order
aaa = a_i(imix+imax-(ii2));
a2 = interp1(ppp,aaa,psis_i(ii1)); 


%%
% interpolation result

plot(a1,psis_i(ii1),'bo',a2,psis_i(ii1),'rx',aaa,ppp,'k.');

%%
a3 = a2 - a1;

% finding position of the island
aa = interp1(a3,a1,w_crit); % a at wcrit
i1 = find(a3<=w_crit,1,'first');
psi_w_crit = interp1(a3,psis_i(ii1),w_crit);
i2 = find(psis_i(ii2) <= psi_w_crit, 1, 'first') + imax - 1; % right side of island

% finding position1 of the island
%i1 = find(a3<=w_crit,1,'first'); % to the left of the island
%psi_w_crit = interp1(a3,psis_i(ii1),w_crit);
%plot(a3,psis_i(ii1),'b.-',w_crit,psi_w_crit,'rx')
%i2 = find(psis_i(ii2) >= psi_w_crit, 1, 'first') + imax - 1; % right side of island

%evaluate quantities at i1 = left side of the island
%aa = interp1(psis_i(ii2),a_i(ii2),psi_w_crit);
rr = interp1(a_i,r_i,aa);
vv = interp1(a_i,vol_i,aa);
nne = interp1(a_i,ne_i,aa);
tte = interp1(a_i,Te_i,aa);

%adding the  values at i2
a_12 = insert_1point(a_i, aa, i2);
r_12 = insert_1point(r_i, rr, i2);
vol_12 = insert_1point(vol_i, vv, i2);
psis_12 = insert_1point(psis_i, psis_i(i1), i2);
ne_12 = insert_1point(ne_i, nne, i2);
Te_12 = insert_1point(Te_i, tte, i2);
lf2 = length(psis_12);

%% Psistar new profile
% [~,psisp] = interpos(r_i,psis_i,tension,bc1,bc2);
psisp = diff(psis_12) ./ diff(a_12);
psis_f(1:i1-1) = psisp(i1)/(2.*a_12(i1)).*(a_12(1:i1-1).^2 - a_12(i1).^2) + psis_12(i1);
psis_f(i1:i2) = psis_12(i2+1);
psis_f((i2+1):lf2) = psis_12((i2+1):lf2);

if saw_params.check
    hs=subplot(222);
    plot(hs,model.rgrid.rhogauss,psis_gauss,'b.-',...
        r_12,psis_f,'rx',...
        r_i(imix), psis_i(imix), 'c.');
    axis([0 1 0 5e-3]);
    ylabel('psis');
    xlabel('r');
    drawnow
end

% calculation of psi profile
% phi1 = pi * B0* r_12.^2 .* rtor_gauss(end) .^ 2; % toroidal flux
phi = pi * B0 * a_12.^2; % toroidal flux
psi_f = psis_f + phi + psi_gauss(1);
psi_f_rhogauss = interp1(r_12, psi_f, rhogauss, 'cubic');
xnew(model.psi.xind) = model.psi.Lamgauss \ psi_f_rhogauss;

if saw_params.check
    hs=subplot(223); hold(hs,'on');
    plot(hs,model.rgrid.rhogauss,1 ./(model.out.Ciotagauss * xnew(model.psi.xind)),'b',... 
    model.rgrid.rhogauss,1./(model.out.Ciotagauss*x(model.psi.xind)),'r');
    xlabel('r');
    ylabel('q');
end

% calculation of density profiles
% [~,~,~,neint] = interpos(vol_i,ne_i,tension,bc1,bc2);
neint = cumtrapz(vol_12, ne_12);
if strcmp(model.ne.method,'state')
    ne_f(1:i2) = neint(i2) ./ vol_12(i2); % reconnected ne profile
else
    ne_f(1:i2) = ne_12(1:i2); % fixed ne profile
end
ne_f(i2+1:lf2) = ne_12(i2+1:lf2);
ne_f_rhogauss = interp1(r_12, ne_f, rhogauss);
if strcmp(model.ne.method,'state')
    xnew(model.ne.xind) = model.ne.LamTegauss \ ne_f_rhogauss;
end

% calculation of temperature profiles
% [~,~,~,peint] = interpos(vol_i,ne_i.*Te_i,tension,bc1,bc2);
peint = cumtrapz(vol_12, ne_12 .* Te_12);
Te_f(1:i2) = peint(i2) ./ (ne_f(i2) * vol_12(i2));
Te_f(i2+1:lf2) = Te_12(i2+1:lf2);
Te_f_rhogauss = interp1(r_12, Te_f, rhogauss);
xnew(model.te.xind) = model.te.LamTegauss \ Te_f_rhogauss;


if saw_params.check
    hs=subplot(224);
    plot(hs,[Te_gauss,Te_f_rhogauss]);
    xlabel('r');
    ylabel('Te');
    drawnow
end

varargout{1} = xnew;
return

function q = insert_1point(q,q1,point1)
q = [q(1:point1),q1,q(point1+1:end)];
return

function q = insert_2points(q,q1,point1,q2,point2)
q = [q(1:point1),q1,q(point1+1:point2),q2,q(point2+1:end)];
return
