Clean up version of running ITER for different cases

IDS/LLXLY_IDS_corsica.m

+function [L,LX,LY] = LLXLY_IDS_corsica(L_LIU, a)
+%% Load data drom IDS
+
+
+% Get geomtrical information from the existing LIUQE
+L = struct();
+L.P = L_LIU.P;
+fieldlist = {'rl';'zl'; 'rw';'zw'; 'Twa'; 'rv';'zv'; 'Tvs'; 'rf'; 'zf';'rm'; 'zm'; 'am'};
+for ii= 1:numel(fieldlist)
+   tmpname = fieldlist{ii};
+   L.G.(tmpname)  = L_LIU.G.(tmpname);
+end
+
+% Set new grid
+L.G.rx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim1;
+L.G.zx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim2;
+
+
+% Recompute the geometrical quantities for the CORSICA grid
+L.G = liug(L.G,L.P);
+L = liuc(L.P, L.G); 
+
+[rrx, zzx] = meshgrid(L.G.rx, L.G.zx);
+
+%%  Get LX structure
+% {Ia, Ip, rBt} available
+% Ff, Bm, Ft recomputed
+
+LX.t = a.pf_active.coil{1}.current.time;
+LX.nt = numel(LX.t);
+
+LX.Ip = zeros(LX.nt,1);
+for ii=1:LX.nt
+    LX.Ip(ii) = a.equilibrium.time_slice{ii}.global_quantities.ip;
+end
+
+% Ia is the current per number of turns
+LX.Ia = zeros(L.G.na, LX.nt);
+for ii= 1:L.G.na
+   LX.Ia(ii,:) = a.pf_active.coil{ii}.current.data/sum(L.G.Twa(:,ii));
+end
+
+LX.rBt = zeros(LX.nt,1);
+LX.rBt(:) = a.equilibrium.vacuum_toroidal_field.r0*a.equilibrium.vacuum_toroidal_field.b0;
+
+
+%% Synthetic data that needs to be recomputed
+LY.Ix = zeros( L.nzx, L.nrx ,LX.nt);
+for tt=1:LX.nt
+    LY.Ix(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor./rrx*L.drx*L.dzx*1e-2; % Wrong definition in IDS Corsica
+   % Ix_Corsica(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor; % Wrong definition of CORSICA
+   % Ix_Corsica(:,:,tt) = -Ix_Corsica(:,:,tt)*abs(LX.Ip(tt)/sum(sum(Ix_Corsica(:,:,tt))));
+end
+
+
+LX.Bm = zeros(L.G.nm, LX.nt);
+for tt= 1:LX.nt
+    LX.Bm(:,tt) = L.G.Bmx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 )  + L.G.Bma*LX.Ia(:,tt);
+end
+
+LX.Ff = zeros(L.G.nf, LX.nt);
+for tt= 1:LX.nt
+    LX.Ff(:,tt) =  L.G.Mfx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 )  + L.G.Mfa*LX.Ia(:,tt);
+end
+
+LX.Is = zeros(0,LX.nt);
+
+LX.shot = 5000;
+
+LY.Fx = zeros(L.nzx,L.nrx ,LX.nt);
+for tt = 1:LX.nt
+   LY.Fx(:,:,tt) =  a.equilibrium.time_slice{tt}.profiles_2d{1}.psi;
+end
+
+% Normalized flux of the R,Z grid
+LY.Opy = zeros(L.nzy, L.nry, LX.nt);
+LY.FNy = zeros(L.nzy, L.nry, LX.nt);
+for tt=1:LX.nt
+[LY.rA(tt),LY.zA(tt),LY.FA(tt),dr2FA,dz2FA,drzFA,rX,zX,~,~,~,~,rl,zl,...
+    LY.rB(tt),LY.zB(tt),LY.FB(tt),LY.lB(tt),LY.lX(tt),LY.Opy(:,:,tt)] = meqpdom(LY.Fx(:,:,tt),LX.Ip(tt),L.P.isaddl,L);
+LY.FNy(:,:,tt) =  (LY.Fx(2:end-1,2:end-1,tt) - LY.FA(tt))/(LY.FB(tt) - LY.FA(tt));
+end
+
+
+
+%% Recompute the plasma current as a check. There is a bit of a difference
+psiN = (a.equilibrium.time_slice{tt}.profiles_1d.psi - a.equilibrium.time_slice{tt}.global_quantities.psi_axis)/(a.equilibrium.time_slice{tt}.global_quantities.psi_boundary - a.equilibrium.time_slice{tt}.global_quantities.psi_axis);
+for tt=1:LX.nt
+pprime = a.equilibrium.time_slice{tt}.profiles_1d.dpressure_dpsi;
+p = spline(psiN,pprime);
+pprimey = ppval(p, LY.FNy(:,:,tt));
+pprimey = pprimey.*LY.Opy(:,:,tt);
+
+ttprime = a.equilibrium.time_slice{tt}.profiles_1d.f_df_dpsi;
+T2d2 = cumtrapz(a.equilibrium.time_slice{tt}.profiles_1d.psi,ttprime);
+T2d2 = T2d2 - T2d2(end);
+
+p = spline(psiN,ttprime);
+ttprimey = ppval(p, LY.FNy(:,:,tt));
+ttprimey = ttprimey.*LY.Opy(:,:,tt);
+
+p = spline(psiN,T2d2);
+Ty = ppval(p, LY.FNy(:,:,tt));
+Ty = Ty.*LY.Opy(:,:,tt);
+
+L.rymap = repmat(L.ry',numel(L.zy),1); % Matriz map ry
+LX.Ft(tt) = sum(sum( Ty/LX.rBt(tt)*L.drx*L.dzx./L.rymap ));
+
+Iy = 2*pi*(L.rymap.*pprimey+1./L.rymap.*ttprimey/mu0)*L.drx*L.dzx;
+
+LY.Ip(tt) = sum(sum(Iy));
+LY.FA(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_axis;
+LY.FB(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_boundary;
+LY.Wk(tt) = a.equilibrium.time_slice{tt}.global_quantities.w_mhd;
+end
+LY.t = LX.t;
+
+
+% Put all vector in raw form
+flist = fieldnames(LX);
+nflist = numel(flist);
+for k = 1:nflist
+    LX.(flist{k}) = v2rw(LX.(flist{k}));
+end
+
+% Put all vector in raw form
+flist = fieldnames(LY);
+nflist = numel(flist);
+for k = 1:nflist
+    LY.(flist{k}) = v2rw(LY.(flist{k}));
+end
+
+
+end
+
+
+% Put all vectors in raw form
+function vout = v2rw(vin)
+if numel(size(vin))==2 && sum(ismember(1,size(vin)))==1 && ~ischar(vin)
+   vout = reshape(vin, numel(vin),1);
+else
+    vout = vin;
+end
+end