imake_equil_ideal

 LUKE - Create an ideal magnetic equilibrium structure for C3PO/LUKE

   Create an ideal magnetic equilibrium structure for C3PO/LUKE. Circular
   and concentric magnetic flux surfaces. 

   by Y. Peysson (CEA/DSM/IRFM,yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM,joan.decker@cea.fr)

0001 function [equil_data,shotnum,shotime] = imake_equil_ideal(e_opt)
0002 %
0003 % LUKE - Create an ideal magnetic equilibrium structure for C3PO/LUKE
0004 %
0005 %   Create an ideal magnetic equilibrium structure for C3PO/LUKE. Circular
0006 %   and concentric magnetic flux surfaces.
0007 %
0008 %   by Y. Peysson (CEA/DSM/IRFM,yves.peysson@cea.fr) and J. Decker (CEA/DSM/IRFM,joan.decker@cea.fr)
0009 %
0010 if nargin < 1,
0011     e_opt = NaN;
0012 end
0013 %
0014 shotnum = '';
0015 shotime = '';
0016 %
0017 if ~isnan(e_opt),% magnetic equilibrium data
0018     %
0019     disp(' ');
0020     info_dke_yp(1,'Manual equilibrium construction : please provide the following data : ')
0021     disp(' ');
0022     %
0023     shotnum = input_dke_yp('Shot number','');
0024     shotime = input_dke_yp('Shot time',[],'','',[1,1]);
0025     %
0026     shotime = num2str(shotime,'%6.4f');
0027     %
0028     R_sep = input_dke_yp('Geometrical major radius of the separatrix (m) (R_sep = Inf for the cylindrical limit)',[],[0;Inf],'',[1,1]);
0029     Z_sep = input_dke_yp('Geometrical vertical shift of the separatrix (m)',0,'','',[1,1]);
0030     a_sep = input_dke_yp('Geometrical minor radius of the separatrix (m)',[],[0;R_sep],'',[1,1]);
0031     %
0032     if e_opt == 1,
0033         %
0034         disp(' ');
0035         info_dke_yp(1,'Parametric description of the separatrix');
0036         disp(' - R_sep = (R_sep_max + R_sep_min)/2, for the guess')
0037         disp(' - a_sep = (R_sep_max - R_sep_min)/2, for the guess.')
0038         disp(' - Z_sep = Z(R_sep_max), for the guess.')
0039         disp(' - elongation_sep = (Z_sep_max - Z_sep_min)/2 (always positive) (Upper triangularity + Lower triangularity in minor radius unit (m) is approximately elongation_sep*a_sep)')
0040         disp(' - upper_triangularity_sep = (R_sep - R(Z_sep_max))/2 (always positive). Upper triangularity in minor radius unit (m) = upper_triangularity_sep*a_sep')
0041         disp(' - lower_triangularity_sep = (R_sep - R(Z_sep_min))/2 (always positive). Lower triangularity in minor radius unit (m) = lower_triangularity_sep*a_sep')
0042         disp('----------------------------------------------------------------------')
0043         disp(' If no upper X-point exists, upper separatrix angles at LFS (R,X) and HFS (-R,X) sides have to be set to 0.')
0044         disp(' If no lower X-point exists, lower separatrix angles at LFS (R,X) and HFS (-R,X) sides have to be set to 0.')
0045         disp(' If an upper X-point exists, upper separatrix angles at LFS (R,X) + upper separatrix angles at LFS (-R,X) is close to 90 degrees')
0046         disp(' If an lower X-point exists, lower separatrix angles at LFS (R,X) + lower separatrix angles at LFS (-R,X) is close to 90 degrees')
0047         disp('----------------------------------------------------------------------')
0048         info_dke_yp(2,'Default values for xpoint correspond to a plasma with a circular poloidal cross-section.');
0049         disp('----------------------------------------------------------------------')
0050         disp(' ');
0051         %
0052         xpoint(1) = input_dke_yp('Upper altitude X point in minor radius unit',1,[0;Inf],'',[1,1]);
0053         xpoint(2) = input_dke_yp('Upper triangularity in minor radius unit (m)',0,[],'',[1,1]);
0054         xpoint(3) = input_dke_yp('Upper separatrix angle (R,X)  (LFS, degrees)',0,[],'',[1,1]);
0055         xpoint(4) = input_dke_yp('Upper separatrix angle (-R,X) (HFS, degrees)',0,[],'',[1,1]);
0056         xpoint(5) = input_dke_yp('Lower altitude X point in minor radius unit (always positive)',1,[0;Inf],'',[1,1]);
0057         xpoint(6) = input_dke_yp('Lower triangularity in minor radius unit',0,[],'',[1,1]);
0058         xpoint(7) = input_dke_yp('Lower separatrix angle (R,X)  (LFS, degrees)',0,[],'',[1,1]);
0059         xpoint(8) = input_dke_yp('Lower separatrix angle (-R,X)  (HFS, degrees)',0,[],'',[1,1]); 
0060         %
0061         equil_data.xpoint = xpoint;
0062     end    
0063     %
0064     Ip = input_dke_yp('Signed value of the plasma current (MA) (positive if clockwise in top view)',[],'','',[1,1]);
0065     Bt = input_dke_yp('Signed value of the toroidal magnetic field on axis (T) (positive if clockwise in top view)',[],'','',[1,1]);
0066     %
0067     qopt = input_dke_yp(['Safety factor profile : ',...
0068         '\n(0) constant (uniform current density, psi=rho^2)',...
0069         '\n(1) profile ofthe form : q = q0 + (qe-q0)*rho^eq'],...
0070         0,{0,1},'Invalid value');
0071     %
0072     if qopt == 1,
0073         %
0074         if isinf(R_sep),
0075             qmin_Rpap = input_dke_yp('Safety factor q0 at plasma center, multiplied by R_sep/a_sep',[],[0;Inf],'',[1,1]);
0076         else
0077             qmin_Rpap = R_sep/a_sep*input_dke_yp('Safety factor q0 at plasma center',[],[0;Inf],'',[1,1]);
0078         end
0079         disp('qe is calculated from Ampere''s law at the plasma edge with a circular plasma cross-section.')
0080         eq = input_dke_yp('Exponent for q radial profile',[],[],'',[1,1]);
0081         %
0082     else
0083         qmin_Rpap = NaN;
0084         eq = NaN;    
0085     end
0086     %
0087     equil_data.Rp = R_sep;
0088     equil_data.Zp = Z_sep;
0089     equil_data.ap = a_sep;
0090     equil_data.Ip = Ip;
0091     equil_data.Bt = Bt;
0092     equil_data.qopt = qopt;
0093     equil_data.qmin_Rpap = qmin_Rpap;
0094     equil_data.eq = eq;
0095     %
0096 else % equilibrium profiles data
0097     %
0098     Zimp = input_dke_yp('Impurity charge (6 for carbon)',6,[],'',[1,1]);
0099     mimp = input_dke_yp('Impurity mass (uma)',12,[],'',[1,1]);
0100     fi = input_dke_yp('Hydrogen isotopic fraction (H/D/T) [1,3]',[0,0.5,0.5],[],'',[1,3]);
0101     %
0102     Te0 = input_dke_yp('Core electron temperature (keV)',[],[],'',[1,1]);
0103     Tea = input_dke_yp('Edge electron temperature (keV)',[],[],'',[1,1]);
0104     eTe1 = input_dke_yp('Exponent #1 for Te profile (Te(r) = (Te0-Tea)*(1-(r/a)^eTe1)^eTe2 + Tea)',[],[],'',[1,1]);
0105     eTe2 = input_dke_yp('Exponent #2 for Te profile (Te(r) = (Te0-Tea)*(1-(r/a)^eTe1)^eTe2 + Tea)',[],[],'',[1,1]);
0106     %
0107     Ti0 = input_dke_yp('Core ion temperature (keV)',[],[],'',[1,1]);
0108     Tia = input_dke_yp('Edge ion temperature (keV)',[],[],'',[1,1]);
0109     eTi1 = input_dke_yp('Exponent #1 for Ti profile (Ti(r) = (Ti0-Tia)*(1-(r/a)^eTi1)^eTi2 + Tia)',[],[],'',[1,1]);
0110     eTi2 = input_dke_yp('Exponent #2 for Ti profile (Ti(r) = (Ti0-Tia)*(1-(r/a)^eTi1)^eTi2 + Tia)',[],[],'',[1,1]);
0111     %
0112     ne0 = input_dke_yp('Core electron density (m-3)',[],[],'',[1,1]);
0113     nea = input_dke_yp('Edge electron density (m-3)',[],[],'',[1,1]);
0114     ene1 = input_dke_yp('Exponent #1 for ne profile (ne(r) = (ne0-nea)*(1-(r/a)^ene1)^ene2 + nea)',[],[],'',[1,1]);
0115     ene2 = input_dke_yp('Exponent #2 for ne profile (ne(r) = (ne0-nea)*(1-(r/a)^ene1)^ene2 + nea)',[],[],'',[1,1]);
0116     %
0117     Zeff0 = input_dke_yp('Core effective charge Zeff (a.u.)',1,[1;Zimp],'',[1,1]);
0118     Zeffa = input_dke_yp('Edge effective charge Zeff (a.u.)',1,[1;Zimp],'',[1,1]);
0119     Zeffopt = input_dke_yp(['Effective charge profile :\n',...
0120         '(0): Zeff(r) = (Zeff0-Zeffa)*(1-(r/a)^2)^eZeff + Zeffa,\n',...
0121         '(1): 0.5*(Zeffa - Zeff0)*tanh(eZeff*(r - xZeff)) + 0.5*(Zeffa + Zeff0)\n'],...
0122         0,{0,1},'Invalid value',[1,1]);
0123     eZeff = input_dke_yp('Exponent for the effective charge profile',[],[],'',[1,1]);
0124     %
0125     if Zeffopt == 1,
0126         xZeff = input_dke_yp('Radial inflexion point for the effective charge profile',[],[0;1],'',[1,1]);
0127     else
0128         xZeff = NaN;
0129     end
0130     %
0131     equil_data.Zi = [1,1,1,Zimp];
0132     equil_data.mi = [1,2,3,mimp];
0133     equil_data.fi = fi;
0134     equil_data.Te0 = Te0;
0135     equil_data.Tea = Tea;
0136     equil_data.eTe1 = eTe1;
0137     equil_data.eTe2 = eTe2;
0138     equil_data.ne0 = ne0;
0139     equil_data.nea = nea;
0140     equil_data.ene1 = ene1;
0141     equil_data.ene2 = ene2;
0142     equil_data.Ti0 = Ti0;
0143     equil_data.Tia = Tia;
0144     equil_data.eTi1 = eTi1;
0145     equil_data.eTi2 = eTi2;
0146     equil_data.Zeff0 = Zeff0;
0147     equil_data.Zeffa = Zeffa;
0148     equil_data.eZeff = eZeff;
0149     equil_data.xZeff = xZeff;
0150     %
0151 end
0152 %
0153