main_C3PO_jd

C3PO - Set ray initial characteristics and launches the C3PO ray-tracing code

 This function calculates the ray initial characteristics and launches the
 C3PO ray-tracing code

 by J. Decker (CEA-IRFM) <joan.decker@cea.fr> and   Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr>

0001 function [wave] = main_C3PO_jd(dkepath,id_wave,equil,equil_fit,rayinit,waveparam,fitparam,rayparam,C3POdisplay,C3POparam,f0struct,fluct,flag)
0002 %C3PO - Set ray initial characteristics and launches the C3PO ray-tracing code
0003 %
0004 % This function calculates the ray initial characteristics and launches the
0005 % C3PO ray-tracing code
0006 %
0007 % by J. Decker (CEA-IRFM) <joan.decker@cea.fr> and   Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr>
0008 %
0009 if nargin < 10,
0010     C3POparam = [];
0011 end
0012 %
0013 if nargin < 11,
0014     f0struct = []; 
0015 end
0016 %
0017 if nargin < 12,
0018     fluct = [];
0019 end
0020 %
0021 if nargin < 13,
0022     flag = 0;%numerical equilibrium used by default
0023 end
0024 %
0025 if ~isempty(fluct) && ~isstruct(fluct),
0026     fluct = [];
0027     disp('WARNING: <fluct> should be a structure or empty in main_C3PO_jd.m. It has been set to [].')
0028 end
0029 %
0030 if ~isempty(f0struct) && ~isstruct(f0struct),
0031     f0struct = [];
0032     disp('WARNING: <f0struct> should be a structure or empty in main_C3PO_jd.m. It has been set to [].')
0033 end
0034 %
0035 if ~isfield(C3POdisplay,'equilibrium'),
0036     C3POdisplay.equilibrium = 0;
0037 end
0038 %
0039 if ~isfield(C3POdisplay,'fluctuations'),
0040     C3POdisplay.fluctuations = 0;
0041 end
0042 %
0043 % Build or load vectorized magnetic equilibrium structure
0044 %
0045 if ~isstruct(equil),
0046     error('The ''equil'' argument in main_C3PO_jd.m must be a structure.') 
0047 end
0048 %
0049 % SO FAR equil IS NEEDED FOR LOOP_MAIN (VESSEL REFLECTIONS).
0050 % Uncomment following lines when only equil_fit will be needed
0051 %
0052 % if isfield(equil,'psin_fit'),%vectorized equilibrium is passed as equil structure
0053 %     equil_fit = equil;
0054 % elseif ~isempty(fitparam),
0055 % ----
0056 if isempty(equil_fit) | length(fieldnames(equil_fit.Bx_fit)) == 7,%no second derivatives for coordinate changes in the plasma core
0057     info_dke_yp(2,['Vectorized form of the magnetic equilibrium ',equil.id,' is being calculated...'],C3POdisplay.equilibrium);
0058     equil_fit = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);
0059     info_dke_yp(2,'    ...done',C3POdisplay.equilibrium);
0060 % ----
0061 % else
0062 %     error('When fitparam structure is empty, the ''equil'' argument in main_C3PO_jd.m must be a structure that corresponds to the interpolated magnetic equilibrium.')
0063 end
0064 %
0065 fluct_fit = [];
0066 if isempty(fluct),
0067     info_dke_yp(2,'No equilibrium fluctuation effects included.',C3POdisplay.fluctuations);
0068 elseif rayparam.tensortype > 0
0069     disp('WARNING: Plasma fluctuations model inconsistent with the plasma dielectric tensor type, fluctuation effects discarded.');
0070 elseif ~isstruct(fluct),
0071     disp('WARNING: The ''fluct'' argument in main_C3PO_jd.m must be a structure, fluctuation effects discarded.') 
0072 elseif isfield(fluct,'xL_perp_fit'), 
0073     fluct_fit = fluct;%reload existing fluctuation structure
0074 elseif isfield(fluct,'xL_perp'), 
0075     if isfield(fluct,'fitparam'),
0076         info_dke_yp(2,['Vectorized form of the plasma fluctuations ',equil.id,'_',fluct.id,' being calculated...'],C3POdisplay.fluctuations);
0077         fluct_fit = fitfluct_yp(fluct,fitparam.mode_equil,fluct.fitparam.method,fluct.fitparam.ngridresample,fluct.fitparam.nharm);
0078         info_dke_yp(2,'    ...done',C3POdisplay.fluctuations);
0079         if C3POdisplay.fluctuations,
0080             testfitfluct_yp(equil_fit,fluct,fluct_fit);
0081         end
0082     else
0083         disp('WARNING: empty field ''fitparam'' in structure ''fluct'', fluctuation effects discarded.')
0084     end
0085 elseif isfield(fluct,'npar0_lh'), 
0086     %to validate that the field exists in the 'fluct' structure
0087 else
0088     disp('WARNING: The ''fluct'' structure is not recognized, fluctuation effects discarded.')
0089 end
0090 %
0091 % wave structure initialization
0092 %
0093 Nrays = length(rayinit.yrho0);
0094 %
0095 omega_rf = rayinit.omega_rf;
0096 %
0097 wave.id = id_wave;
0098 wave.model = 'RT';
0099 wave.equil_id = equil.id;
0100 wave.omega_rf = omega_rf;
0101 wave.opt_rf = waveparam.opt_rf;
0102 if isfield(waveparam,'n_rf_list'),
0103     wave.n_rf_list = waveparam.n_rf_list;
0104 end
0105 %
0106 wave.opt_disp = 0;  % option for solving the wave equation, which gives (ray.sepol_pmz, ray.sphi, ray.sNperp)
0107 %
0108 wave.dsmin = waveparam.dsmin;
0109 %
0110 % Backward compatibility
0111 %
0112 if ~isfield(C3POdisplay,'text'),
0113     C3POdisplay.text = 1;
0114 end
0115 %
0116 if ~isfield(C3POdisplay,'mdce'),
0117     C3POdisplay.mdce = 0;
0118 end
0119 %
0120 if ~isfield(waveparam,'nd'),
0121     waveparam.nd = 1;
0122 end
0123 if ~isfield(waveparam,'nchi'),
0124     waveparam.nchi = 1;
0125 end
0126 if ~isfield(waveparam,'ns'),
0127     waveparam.ns = 1;
0128 end
0129 if isfield(waveparam,'a_sdNpar'),
0130     wave.a_sdNpar = waveparam.a_sdNpar;
0131 end
0132 %
0133 if ~isfield(rayparam,'dS'),% !!!! old rayparam, modify calling function
0134     disp('WARNING: You should update old rayparam')
0135     rayparam.dS = rayparam.drho;
0136     %
0137     rayparam.rel_opt = 1;%option for (1) relativistic or (0) non-relativistic calculations
0138     rayparam.nperp = 1000;%number of points in pperp integration for damping calculations
0139     rayparam.pperpmax = 10;%maximum value of pperp in damping calculations
0140     rayparam.tau_lim = 20;%value of optical depth beyond which the wave is considered absorbed
0141     %
0142     rayparam = rmfield(rayparam,'drho');
0143 end
0144 %
0145 if ~isfield(rayparam,'kextra'),% !!!! old rayparam, modify calling function
0146     rayparam.kextra = 0;
0147 end
0148 %
0149 if ~isfield(rayparam,'rhomin'),% minimum distance of the ray to the magnetic axis
0150     rayparam.rhomin = 0;
0151 end
0152 %
0153 if ~isfield(rayparam,'metricmode'),% 1 : Only curvilinear metric; 2 : Both curvilinear and cartesian metric
0154     rayparam.metricmode = 1;
0155 end
0156 %
0157 if ~isfield(rayparam,'rhoswitch'),% Radial position at which the switch occurs when metricmode = 2
0158     rayparam.rhoswitch = 0.1;
0159 end
0160 %
0161 if ~isfield(rayparam,'deltaswitch'),% Delta  switch between (rho,theta) <-> (X,Y) to avoid oscillations between modes
0162     rayparam.deltaswitch = 0.01;
0163 end
0164 %
0165 if ~isfield(rayparam,'substitution_method'),% %Method of substitution : 1->dichotomy1D, 2->dichotomy2D, 3->newton_numeric1D, 4->newton_numeric2D, 5->newton_analytic1D, 6->newton_analytic2D, 7->newton_mix2D
0166     rayparam.substitution_method =7;
0167 end
0168 %
0169 if flag == 2 && rayparam.metricmode == 2;
0170     flag = 2;
0171     rayparam.metricmode = 1;
0172     disp('WARNING: analytic tests must use (rho,theta) metric: rayparam.metricmode = 1 enforced.');
0173 end
0174 %
0175 %if ~isfield(rayinit,'yn0'),% !!!! old rayparam, modify calling function
0176 %    rayinit.yn0 = NaN*rayinit.ym0;
0177 %end
0178 %
0179 %if ~isfield(rayinit,'ykz0'),% !!!! old rayparam, modify calling function
0180 %    rayinit.ykz0 = NaN*rayinit.ym0;
0181 %end
0182 %
0183 if ~isfield(rayinit,'w0'),% waist information
0184     rayinit.w0 = NaN;
0185 end
0186 if ~isfield(rayinit,'z0'),% waist position information
0187     rayinit.z0 = NaN;
0188 end
0189 %
0190 waveparam.cut_fac = rayparam.tau_lim/10;
0191 %
0192 wave.cut_fac = waveparam.cut_fac;
0193 wave.kextra = rayparam.kextra;
0194 %
0195 arg_range = cell(1,Nrays);
0196 for iray = 1:Nrays,
0197     arg_range{iray} = iray;
0198 end
0199 %
0200 if isempty(dkepath),
0201     [dkepath] = load_structures_yp('dkepath','','');
0202 end
0203 %
0204 if isempty(C3POparam),
0205     C3POparam.clustermode.main_C3PO_jd.scheduler.mode = 0;%sequential mode
0206 end
0207 %
0208 dkecluster = clustermode_luke(C3POparam.clustermode,'main_C3PO_jd',dkepath);%MatLab distributed computing environment
0209 %
0210 [flag_mdce_main_C3PO_jd,rayss] = mdce_run(@loop_main_C3PO_jd,{iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,waveparam,C3POdisplay,f0struct,fluct_fit},1,arg_range,dkecluster);
0211 %
0212 iray_wave = 1;
0213 for iray = 1:Nrays,
0214     %
0215     rays = rayss{iray};
0216     Niray = length(rays);%ray can be duplicated if multiple reflections
0217     %
0218     if Niray == 1 && isempty(rays{1}.srho),
0219         disp(['kpsin has a non nul imaginary part for ray # ',int2str(iray),' !']);  
0220         info_dke_yp(2,['WARNING: Ray #',int2str(iray),' is not put in the wave structure.']);
0221     else
0222         wave.rays(iray_wave:iray_wave + Niray - 1) = rays;
0223         iray_wave = iray_wave + Niray;
0224         info_dke_yp(2,['Ray #',int2str(iray),' is processed']);    
0225     end
0226 end    
0227 clear rayss rays
0228 %
0229 wave = make_beam_wave_jd(wave,equil,C3POdisplay,waveparam.nd,waveparam.nchi,rayinit.w0,rayinit.z0,waveparam.ns);
0230 %
0231 if length(rayinit.yP0_2piRp) > length(wave.rays),
0232     ptot_coupled = 0;
0233     %
0234     for ii = 1:length(wave.rays)
0235         ptot_coupled = ptot_coupled + wave.rays{ii}.P0_2piRp*2*pi*equil.Rp;
0236     end
0237     %
0238     disp('WARNING: Some rays have been removed since they do not penetrate the plasma (kpsin has a non nul imaginary part)');
0239     disp(['WARNING: The launched rf power is : ',num2str(sum(rayinit.yP0_2piRp)*2*pi*equil.Rp/1e6),' (MW), while the rf power coupled to the plasma is : ',num2str(ptot_coupled/1e6),' (MW).']);
0240 end
0241 %
0242 wave.rayinit = rayinit;
0243 wave.waveparam = waveparam;
0244 wave.fitparam = fitparam;
0245 wave.rayparam = rayparam;
0246 wave.C3POparam = C3POparam;
0247 wave.display = C3POdisplay;
0248 wave.flag_mdce_main_C3PO_jd = flag_mdce_main_C3PO_jd;
0249 %
0250