function Parameters = fgetk_environment_parameters(testid, shot, time, reload,doplot)
% Get parameters to be used by fgetk_environment.m
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Load or Retrieve/Save data
tok = 'tcv';
idoublet = false; % default
switch testid
  case 'nl_cde'
    lin = false;
    agcon = {'cde0Dss','bp','qA'};
  case 'lin_cde'
    lin = true;
    agcon = {'cde0Dss','bp','qA'};r
  case 'nl_fixIp'
    lin = false;
    agcon = {'Ip','bp','qA'};
  case 'lin_fixIp'
    lin = true;
    agcon = {'Ip','bp','qA'};
  case 'droplet_fixIp'
    lin = 'false';
    idoublet = true;
    agcon = {{'Ip'},{'Ip'},{'Ip'}};
  case 'droplet_cde'
    lin = 'false';
    idoublet = true; 
    agcon = {{'cde0Dss'},{'cde0Dss'},{'Ip'}};
  case 'doublet_cde'
    lin = 'false';
    idoublet = true; 
    agcon = {{'cde0Dss'},{'cde0Dss'},{'Ip'}};
  otherwise
    error('unknown testid');
end

if idoublet
  shot = 67151; t1 = 0.03; % overwrite with specific droplet prepared shot
  time = t1+(0:1e-4:1e-3);
end

debug = 0;
iterq = 20;


fpath   = fileparts(mfilename('fullpath'));
fname   = fullfile(fpath,sprintf('fge_TCV%d_t%1.3f_%s.mat',shot,time(1),testid));

if reload && ~isempty(which(fname))
  fprintf('load data from %s\n',fname);
  Parameters = load(fname);
else
  
  if ~idoublet
    [L,LX] = fge(tok,shot,time,'debug',debug,...
      'lin',lin,...
      'agcon',agcon,...
      'izgrid',true,...
      'infct',@qintmex,...
      'iterq',iterq);
    LX = fgex(tok,time,L,LX);
    [L,LX]   = fgel(L,LX); % linearization + initial condition
    
  else
    %% special loader using FBT for doublets
    % Compute equilibrium coil currents using fbt
    [Lfbt,LXfbt] = fbt(tok,shot,time(1),...
      'selx','X','izgrid',true,...
      'issym',0,... % allow asymmetry
      'bfct',@bfabmex,'bfp',[0,1],... % one basis function, zero pressure
      'agfitfct',@meqfit3,'agfitp',[true,false,false],... % constrain only Ip
      'tol',1e-9,'niter',100,'debug',debug,...
      'iterq',iterq);
    
    %% Customizations for droplet positions
     k = ~isnan(LXfbt.gpz); % valid points
     np = sum(k)/2; % number of points per droplet

     if contains(testid,'droplet')
       % play with parameters for different different droplet size/position
       r0 = 0.83; z0 = 0.4; a  = 0.21;  % play with these
       kt = 1.0;  kb = kt; dt = 0.2; db = dt; lt = 0.0; lb = lt; % play with these
       rt = r0+0.05; rb = r0-0.02; zt = z0; zb = -z0-0.05;
       [ra1,za1] = bouD(np,rt,zt,a,kt,kt,dt,dt,lt,lt,0,0); % top droplet
       [ra2,za2] = bouD(np,rb,zb,a,kb,kb,db,db,lb,lb,0,0); % bot droplet
       
       LXfbt.IpD(1:2) = [-190e3,-170e3]; % assign custom/different Ip
       LXfbt.gpr(k) = [ra1;ra2]; % assign new r positions of control points
       LXfbt.gpz(k) = [za1;za2]; % assign new z positions of control points
     elseif contains(testid,'doublet')
       % doublet !
       r0 = 0.83; z0 = 0.25; a  = 0.21;  % play with these
       kt = 1.2;  kb = kt; dt = 0.2; db = dt; lt = 0.0; lb = lt; % play with these
       rt = r0; rb = r0; zt = z0; zb = -z0-0.05;
       [ra1,za1] = bouD(np,rt,zt,a,kt,kt,dt,dt,lt,lt,0,0); % top droplet
       [ra2,za2] = bouD(np,rb,zb,a,kb,kb,db,db,lb,lb,0,0); % bot droplet
       
       LXfbt.IpD(1:2) = [-100e3,-100e3]; % assign custom/different Ip
       LXfbt.gpr(k) = [ra1;ra2]; % assign new r positions of control points
       LXfbt.gpz(k) = [za1;za2]; % assign new z positions of control points
     else, error('this should not happen')
     end
     % meqgplot(Lfbt.G); hold on; plot(ra,za,'o'); % check points
     
    %% compute initial equilibrium and coil currents
    LYfbt = fbtt(Lfbt,LXfbt,'debug',0); 
    assert(~isempty(LYfbt) && LYfbt.isconverged,'FBT did not converge');
    
    if doplot
      meqplott(Lfbt,LYfbt); hold on;
      %plot([ra1;ra2],[za1;za2],'o'); % plot points
      title('Initial FBT equilibrium'); drawnow; 
    end
    
    %% Setup FGE
    L   = fge(tok,shot,time,...
      'idoublet',idoublet,...
      'agcon',agcon,...
      'selu','v','nu',256,...
      'usepreconditioner',1,...
      'bfct',@bfabmex,'bfp',[0,1],...
      'mkryl',20,...
      'issym',false,...
      'izgrid',true,...
      'debug',1,...
      'infct',@qintmex,...
      'iterq',iterq);
    %%
    
    LXin = meqxconvert(Lfbt,LYfbt,L); % initial LX
    [L,LX]   = fgel(L,LXin); % linearization
    LX = fgex('tcv',time,L,LX); % put in fge form

    LX.Va(1,2:10) = LX.Va(1,2:10)+100; % perturb Va a bit so it moves.
    
    %% Customizations for droplet resistivity 
    LX.signeo(1,:) = LX.signeo(1,:)*0.8; % lower conductivity in first droplet
    
    if doplot, clf; meqplotfancy(L,LXin); title('Initial equilibrium'); drawnow; end
  end
end

%% Add points on which we'll interpolate flux, fields - to be used for shape rewards
% points on which intepolated flux, fields are returned
r0 = 0.86; a = min(r0-min(L.G.rl),max(L.G.rl)-r0); % make sure we limit on a wall
[rn,zn] = bouD(16,r0,0,a,1.5,1.5,0.2,0.2,0,0,0,0); 
% uncomment this to plot the points
meqgplot(L.G); hold on; plot(rn,zn,'x')
L.P.rn = rn; L.P.zn = zn; % store in L.P
%%

Parameters.L  = L;
Parameters.LX = LX;
Parameters.MeasParams = LY2Ain(L);

% also save the Khybrid output names
% we use this in various places in the environments
% this replicates functioniality from rzip/Tools/get_mmat.m;
mdsopen('pcs',shot);
mmat.data = mdsdata('\pcs::phys_mat_m');
Parameters.Khybrid.OutputName = mdsdata('\pcs::phys_mat_m_outputs');
fprintf('saved Parameters to %s\n',fname);

save(fname,'-struct','Parameters');
end
