make_equil_metis2luke

 creates LUKE equilibrium structure from METIS output

 MAGNETIC EQUILIBRIUM

0001 function [equil,profil0d] = make_equil_metis2luke(equil_id,profil0d_in,zerod,exp0d,cons,geo,option,equilparam)
0002 %
0003 % creates LUKE equilibrium structure from METIS output
0004 %
0005 % MAGNETIC EQUILIBRIUM
0006 %
0007 if equilparam.helena >= 1,
0008     %
0009     % ***** 2D equilibrium reconstructed by Helena *****
0010     %
0011     post.profil0d = profil0d_in;
0012     post.zerod = zerod;
0013     post.z0dinput.exp0d = exp0d;
0014     post.z0dinput.cons = cons;
0015     post.z0dinput.geo = geo;
0016     post.z0dinput.option = option;
0017     %
0018     if equilparam.helena == 1,% use helena to recalculate equilibrium
0019         %
0020         equi = z0dhelena(post,0);
0021         equil = make_equil_magnetic_cronos(equi);
0022         %
0023         xli = equi.a/equi.a(end);
0024         Raxe = equi.raxe;
0025         kx = equi.e;
0026         dx = equi.d;
0027         %
0028     else
0029         %
0030         [filename,pathname] = uigetfile('EQUIL_*.mat','Choose equilibrium file');
0031         %
0032         load([pathname,filename],'equil');
0033         %
0034         if ~isfield(equil,'pli'),
0035             disp('warning : pli data not available, using circular low beta approximation')
0036             xli = (equil.ptx(:,1) - equil.ptx(:,(end+1)/2))/(equil.ptx(end,1) - equil.ptx(end,(end+1)/2));
0037         else
0038             xli = equil.pli/equil.pli(end);
0039         end
0040         %
0041         Raxe = NaN;
0042         kx = NaN;
0043         dx = NaN;
0044         %
0045     end
0046     %
0047     profil0d.xli = xli;
0048     profil0d.tep = interp1(profil0d_in.xli,profil0d_in.tep,profil0d.xli,'pchip','extrap');
0049     profil0d.tip = interp1(profil0d_in.xli,profil0d_in.tip,profil0d.xli,'pchip','extrap');
0050     profil0d.nep = interp1(profil0d_in.xli,profil0d_in.nep,profil0d.xli,'pchip','extrap');
0051     profil0d.n1p = interp1(profil0d_in.xli,profil0d_in.n1p,profil0d.xli,'pchip','extrap');
0052     profil0d.nhep = interp1(profil0d_in.xli,profil0d_in.nhep,profil0d.xli,'pchip','extrap');
0053     profil0d.nzp = interp1(profil0d_in.xli,profil0d_in.nzp,profil0d.xli,'pchip','extrap');
0054     %
0055     profil0d.epar = interp1(profil0d_in.xli,profil0d_in.epar,profil0d.xli,'pchip','extrap');
0056     %
0057     % trois lignes suivantes à vérifier
0058     %
0059     profil0d.Raxe = Raxe;
0060     profil0d.kx = kx;
0061     profil0d.dx = dx;
0062     %
0063     % sign management
0064     %
0065     sigmaB = sign(equil.ptBPHI(1,1));
0066     sigmaI = sign(equil.ptBy(end,1));
0067     %
0068     equil.ptBPHI = sigmaB*option.bsigne*equil.ptBPHI;
0069     equil.ptBx = sigmaI*option.bsigne*option.signe*equil.ptBx;
0070     equil.ptBy = sigmaI*option.bsigne*option.signe*equil.ptBy;
0071     equil.psi_apRp = sigmaI*option.bsigne*option.signe*equil.psi_apRp;
0072     %
0073 else
0074     %
0075     % ***** 2D equilibrium from moments by metis2equi1t *****
0076     %
0077     nt = 65;% number of poloidal grid points
0078     %
0079     % see function mapequilibrium for call to metis2equi1t
0080     %
0081     % constante physique (phys)
0082     phys.c           =   2.99792458e8;             % vitesse de la lumiere dans le vide (m/s)  (definition)
0083     phys.h           =   6.62606876e-34;           % constante de Planck (J*s) (+/- 0.0000052e-34)
0084     phys.e           =   1.602176462e-19;          % charge de l'electron (C)   (+/- 0.000000063e-19)
0085     phys.mu0         =   4*pi*1e-7;                % permeabilite du vide (H/m) (definition)
0086     phys.epsi0       =   1./phys.c.^2./phys.mu0;   % permitivite du vide (F/m)  (definition)
0087     phys.g           =   6.673e-11;                % constante de la gravitation (N*m^2/kg^2) (+/- 0.010e-11)
0088     phys.k           =   1.3806503e-23;            % constante de Boltzmann (J/K)  (+/- 0.0000024e-23)
0089     phys.alpha       =   7.297352533e-3 ;          % constante de structure fine (+/- 0.000000027e-3 )
0090     phys.me          =   9.10938188e-31;           % masse au repos de l'electron (kg) (+/- 0.00000079e-31)
0091     phys.mp          =   1.6726485e-27;            % masse au repos du proton (kg)
0092     phys.ua          =   1.66053873e-27;           % 1 unite atomique (kg) (+/- 1.00000013e-27)
0093     phys.avo         =   6.02214199e23;            % nombre d'avogadro (mol^-1) (+/- 0.00000047e23)
0094     phys.sigma       =   5.670400e-8;              % constante de stephan ( W*m^-2*K^-4) (+/- 0.000040e-8)
0095     %
0096     % spread equilibrium data on interpolated grid
0097     %
0098     if equilparam.nbpx > 0,
0099         profil0d = z0dinterpx(profil0d_in,equilparam.nbpx,equilparam.mode);
0100     else
0101         profil0d = profil0d_in;
0102     end
0103     %
0104     prof.x    = profil0d.xli;
0105     prof.kx   = profil0d.kx;     
0106     prof.dx   = profil0d.dx;      
0107     prof.rmx  = profil0d.rmx;     
0108     prof.Raxe = profil0d.Raxe;
0109     prof.psi  = profil0d.psi;
0110     prof.fdia = profil0d.fdia;
0111     prof.jmoy = profil0d.jli;
0112     prof.ptot = profil0d.ptot;
0113     %
0114     if isfield(profil0d,'Rsepa') && isfield(profil0d,'Zsepa') && all(isfinite(profil0d.Rsepa(:))) && all(isfinite(profil0d.Zsepa(:)))
0115         geo.Rsepa = profil0d.Rsepa;
0116         geo.Zsepa = profil0d.Zsepa;
0117     else
0118         geo.Rsepa = [];
0119         geo.Zsepa = []; 
0120         equilparam.expo = 0;
0121     end
0122     %
0123     %[profil,deuxd,moment,scalaire,facteur] = metis2equi1t(profil,geo,phys,nbth,facteur,exposhape,regularisation);
0124     %
0125     % facteur is a dummy argument for now (02/01/2011)
0126     %
0127     [profil,deuxd] = metis2equi1t(prof,geo,phys,nt,1,equilparam.expo,0);
0128     %
0129     % In METIS, deuxd.BPHI is always positive, the
0130     % relative sign of Ip is prescribed by option.signe
0131     %
0132     equil.Rp = deuxd.R(1,1);
0133     equil.Zp = deuxd.Z(1,1);
0134     %
0135     equil.ptx = deuxd.R - equil.Rp;
0136     equil.pty = deuxd.Z - equil.Zp;
0137     equil.ptBx = - option.bsigne*option.signe*deuxd.BR;
0138     equil.ptBy = - option.bsigne*option.signe*deuxd.BZ;
0139     equil.ptBPHI = option.bsigne*deuxd.BPHI;
0140     %
0141     equil.psi_apRp = - option.bsigne*option.signe*(profil.psi - profil.psi(1))*equil.ptx(end,1)/equil.Rp;
0142     equil.theta = linspace(0,2*pi,nt);
0143     %
0144 end
0145 %
0146 if equil.psi_apRp(end)*equil.ptBy(end,1) < 0,
0147     disp('-------> WARNING : equilibrium signs not self-consistent')
0148 end
0149 %
0150 if equil.psi_apRp(end)*equil.ptBPHI(1,1)*option.signe < 0,
0151     disp('-------> WARNING : equilibrium signs not consistent with option.signe')
0152 end
0153 %
0154 if equil.ptBPHI(1,1)*option.bsigne < 0,
0155     disp('-------> WARNING : equilibrium signs not consistent with option.bsigne')
0156 end
0157 %
0158 % EQUILIBRIUM PROFILES
0159 %
0160 % see function ./zerod/z0dacces
0161 %
0162 equil.pTe = profil0d.tep/1e3;
0163 equil.pne = profil0d.nep;
0164 %
0165 zimp = fix(option.zimp);
0166 aimp = fix(zimp*7/3);
0167 zmax = fix(option.zmax);
0168 amax = fix(zmax*7/3);
0169 %
0170 if strcmp(option.mino,'He3'),
0171     ahe = 3;
0172 else
0173     ahe = 4;
0174 end
0175 %
0176 equil.zZi = [1,1,1,2,zimp,zmax];
0177 equil.zmi = [1,2,3,ahe,aimp,amax];
0178 %
0179 fD = zerod.nDm/zerod.n1m;
0180 fT = zerod.nTm/zerod.n1m;
0181 fH = max([0,1-fD-fT]);
0182 %
0183 equil.fi = [fH,fD,fT];
0184 pzni1 = equil.fi.'*profil0d.n1p;
0185 %
0186 pnhe = profil0d.nhep;
0187 pnimp = profil0d.nzp;
0188 pnmax = profil0d.nzp*option.rimp;
0189 %
0190 equil.pzni = [pzni1;pnhe;pnimp;pnmax];
0191 %
0192 equil.pzTi = ones(length(equil.zZi),1)*profil0d.tip/1e3;
0193 %
0194 % calculation of moments from the LUKE equilibrium
0195 %
0196 [equil.pshift,equil.pelong,equil.ptrian,equil.pli] = equilmoments_jd(equil.theta,equil.ptx,equil.pty);
0197 %
0198 equil.ip = zerod.ip;% plasma current in A
0199 equil.iohm = zerod.iohm;% ohmic current in A
0200 equil.iboot = zerod.iboot;% bootstrap current in A
0201 %
0202 % modify after error correction in METIS
0203 %
0204 equil.icd = zerod.ipar - zerod.iohm - equil.iboot;% zerod.icd;% zerod.ini - equil.iboot;
0205 %
0206 if equilparam.comp,
0207     %
0208     if equilparam.helena == 'y',
0209         leg = {'METIS - orig.','HELENA','LUKE'};
0210     else
0211         leg = {'METIS - orig.','METIS - interp.','LUKE'};
0212     end
0213     %
0214     figure(101),clf
0215     plot(profil0d_in.xli,profil0d_in.Raxe - equil.Rp,'b-o',profil0d.xli,profil0d.Raxe - equil.Rp,'m-+',equil.pli,equil.pshift,'r--s')
0216     title('Shafranov shift (m)');legend(leg);
0217     %
0218     figure(102),clf
0219     plot(profil0d_in.xli,profil0d_in.kx,'b-o',profil0d.xli,profil0d.kx,'m-+',equil.pli,equil.pelong,'r--s')
0220     title('Elongation');legend(leg);
0221     %
0222     figure(103),clf
0223     plot(profil0d_in.xli,profil0d_in.dx,'b-o',profil0d.xli,profil0d.dx,'m-+',equil.pli,equil.ptrian,'r--s')
0224     title('Triangularity');legend(leg);
0225     %
0226 end
0227 %
0228 equil.id = equil_id;
0229 %
0230 equil = equilconsistency_yp(equil);
0231 %
0232 
0233 
0234 
0235 
0236 
0237 
0238 
0239 
0240 
0241 
0242 
0243 
0244 
0245 
0246 
0247 
0248 
0249 
0250 
0251 
0252 
0253 
0254 
0255