make_equil_cronos

 This function creates the equil structure used by LUKE from various
 cronos structures.

 by Joan Decker (joan.decker@cea.fr) and Yves Peysson (yves.peysson@cea.fr)

0001 function equil = make_equil_cronos(equil_id,equi,gene,prof,compo,impur)
0002 %
0003 % This function creates the equil structure used by LUKE from various
0004 % cronos structures.
0005 %
0006 % by Joan Decker (joan.decker@cea.fr) and Yves Peysson (yves.peysson@cea.fr)
0007 method = 'spline';
0008 %
0009 % protection pour des flux non monotones
0010 % convention : equi.psiRZ(end) - equi.psiRZ(1) < 0
0011 %
0012 rhon_RZ = double(equi.rhoRZ/equi.rhoRZ(end));
0013 dpsi_RZ = diff(double(equi.psiRZ));
0014 %
0015 dpsi_RZ_min = max(dpsi_RZ(dpsi_RZ < 0));
0016 %
0017 dpsi_RZ(dpsi_RZ > 0) = dpsi_RZ_min;
0018 %
0019 psiRZ = cumsum([0,dpsi_RZ]);
0020 %
0021 equil.id = equil_id;%gene.k is the time slice
0022 %
0023 if any(imag(equi.R(:)) ~= 0) || any(imag(equi.Z(:)) ~= 0),
0024     disp('WARNING : imaginary data from Helena equilibrium, discarded in LUKE');
0025 end
0026 %
0027 equil.Rp = double(real(equi.R(1,1,1)));
0028 equil.Zp = double(real(equi.Z(1,1,1)));
0029 equil.ptx = squeeze(double(real(equi.R))) - equil.Rp;
0030 equil.pty = squeeze(double(real(equi.Z))) - equil.Zp;
0031 equil.ap = equil.ptx(end,1);
0032 equil.psi_apRp = -psiRZ*equil.ap/equil.Rp;%Match already the defintion of ppsi in DKE (psi(helena)/2/pi))
0033 equil.theta = linspace(0,2*pi,size(equi.R,3));
0034 equil.ptBx = squeeze(double(equi.BR));
0035 equil.ptBy = squeeze(double(equi.BZ));
0036 equil.ptBPHI = squeeze(double(equi.BPHI));
0037 %
0038 equil.ptBPHI = -gene.signe.b0*equil.ptBPHI;%Bphi is counter-clockwise in TS (WARNING: in rtluke, Ip > 0 and Bphi > 0 when clockwise from the top view of the tokamak, opposite from CRONOS def)
0039 equil.ptBx = -gene.signe.ip*equil.ptBx;%Ip is counter-clockwise in TS
0040 equil.ptBy = -gene.signe.ip*equil.ptBy;%Ip is counter-clockwise in TS
0041 equil.psi_apRp = -gene.signe.ip*equil.psi_apRp;%Ip is counter-clockwise in TS
0042 %
0043 equil.pTe = interp1(gene.x,prof.te/1000,rhon_RZ,method);%en keV
0044 equil.pne = interp1(gene.x,prof.ne,rhon_RZ,method);% m-3
0045 %
0046 zH = find(compo.a == 1);
0047 zD = find(compo.a == 2 & compo.z == 1);
0048 zT = find(compo.a == 3 & compo.z == 1);
0049 zimp = find(compo.z > 1);
0050 %
0051 % profil de densite electronique (m-3)
0052 %wavestructs
0053 pzni = zeros((3+length(zimp)),length(prof.ne));
0054 pTi = interp1(gene.x,prof.ti/1000,rhon_RZ,method);
0055 pzTi = ones((3+length(zimp)),1)*pTi;%All impurities have the same ion temperature
0056 if ~isempty(zH),
0057     pzni(1,:) = interp1(gene.x,squeeze(impur.impur(1,:,zH)),rhon_RZ,method);
0058 end
0059 if ~isempty(zD),
0060     pzni(2,:) = interp1(gene.x,squeeze(impur.impur(1,:,zD)),rhon_RZ,method);
0061 end
0062 if ~isempty(zT),
0063     pzni(3,:) = interp1(gene.x,squeeze(impur.impur(1,:,zT)),rhon_RZ,method);
0064 end
0065 %
0066 pzni(4:(3+length(zimp)),:) = squeeze(impur.impur(1,:,zimp))';
0067 pzTi(4:(3+length(zimp)),:) = repmat(pTi,length(zimp),1);%en keV
0068 zZi = [1,1,1,compo.z(zimp)];
0069 zmi = [1,2,3,compo.a(zimp)];
0070 %
0071 equil.pzTi = pzTi;
0072 equil.pzni = pzni;
0073 equil.zZi = zZi;
0074 equil.zmi = zmi;
0075 %
0076 if sum(equil.ptBy(:,end)-equil.ptBy(:,1)) ~= 0
0077     disp('equilibrium warning, By(:,theta=0) ~= By(:,theta=pi)')
0078     equil.ptBy(:,end) = equil.ptBy(:,1);
0079 end
0080 if sum(equil.ptBx(:,end)-equil.ptBx(:,1)) ~= 0
0081     disp('equilibrium warning, Bx(:,theta=0) ~= Bx(:,theta=pi)')
0082     equil.ptBx(:,end) = equil.ptBx(:,1);
0083 end
0084 if sum(equil.ptBPHI(:,end)-equil.ptBPHI(:,1)) ~= 0
0085     disp('equilibrium warning, Bphi(:,theta=0) ~= Bphi(:,theta=pi)')
0086     equil.ptBPHI(:,end) = equil.ptBPHI(:,1);
0087 end
0088