loop_main_C3PO_jd

C3PO - Set ray initial characteristics and launches the C3PO ray-tracing code (loop function for distributed computing environment)

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

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

 this function is called for a single ray

0001 function [rays,rhomin] = loop_main_C3PO_jd(iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,waveparam,C3POdisplay,f0struct,fluct_fit)
0002 %C3PO - Set ray initial characteristics and launches the C3PO ray-tracing code (loop function for distributed computing environment)
0003 %
0004 %This function calculates the ray initial characteristics and launches the C3PO ray-tracing code
0005 %
0006 %by Yves Peysson (CEA-IRFM,yves.peysson@cea.fr) and Joan Decker (CEA-IRFM,joan.decker@cea.fr) and
0007 %
0008 % this function is called for a single ray
0009 %
0010 if isfield(waveparam,'rhomin') && isfield(rayinit,'type') && strcmp(rayinit.type,'EC'),
0011    %
0012    [rays,rhomin] = loop_main_C3PO_jd(iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,rmfield(waveparam,'rhomin'),C3POdisplay,f0struct,fluct_fit);
0013    %
0014    if rhomin < waveparam.rhomin,
0015         %
0016         disp(' ');
0017         disp(['Minimum rho = ',num2str(rhomin),' below the threshold ',num2str(waveparam.rhomin),' : adjusting rayinit.ym0 to reach threshold']);
0018         %
0019         optim.tolfun = 1e-1;% tolerance on targets
0020         optim.tolx = 1e-1;% tolerance on inputs
0021         optim.nitmax = 20;%
0022         %
0023         ym0 = fminsearch(@(X) calc_rhomin_jd(X,...
0024             iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,waveparam,C3POdisplay,f0struct,fluct_fit),...
0025             rayinit.ym0,optimset('TolFun',optim.tolfun,'TolX',optim.tolx,'MaxIter',optim.nitmax));
0026         %
0027         disp(' ');
0028         disp(['  --> rayinit.ym0 changed from ',num2str(rayinit.ym0),' to ',num2str(ym0)]);
0029         %
0030         rayinit.ym0 = ym0;
0031         %
0032         [rays,rhomin] = loop_main_C3PO_jd(iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,rmfield(waveparam,'rhomin'),C3POdisplay,f0struct,fluct_fit);
0033         %
0034         disp(' ');
0035         disp(['  --> New minimum rho = ',num2str(rhomin)]);
0036         %
0037    end
0038    %
0039    return
0040    %
0041 end
0042 %
0043 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0044 %
0045 rayinit.yrho0 = rayinit.yrho0(iray);
0046 rayinit.ytheta0 = rayinit.ytheta0(iray);
0047 %
0048 if isfield(rayinit,'yphi0') && ~isinf(equil_fit.Rp),
0049     rayinit.yz0 = rayinit.yphi0*equil_fit.Rp;
0050     rayinit = rmfield(rayinit,'yphi0');
0051 elseif isfield(rayinit,'yphi0') && isinf(equil_fit.Rp),
0052     error('ERROR in loop_main_C3PO_jd.m: inconsistent rayinit input: use yphi0 if Rp is finite only, or yz0 elsewhere.');
0053 end
0054 %
0055 rayinit.yz0 = rayinit.yz0(iray);
0056 %
0057 rayinit.ym0 = rayinit.ym0(iray);
0058 %
0059 if isfield(rayinit,'yn0') && ~isinf(equil_fit.Rp),
0060     rayinit.ykz0 = rayinit.yn0/equil_fit.Rp;%backward compatibility for the toroidal coordinates
0061     rayinit = rmfield(rayinit,'yn0');
0062 elseif isfield(rayinit,'yn0') && isinf(equil_fit.Rp),
0063     error('ERROR in loop_main_C3PO_jd.m: inconsistent rayinit input: use yn0 if Rp is finite only, or ykz0 elsewhere.');
0064 end
0065 %
0066 rayinit.ykz0 = rayinit.ykz0(iray);
0067 %
0068 rayinit.yNpar0 = rayinit.yNpar0(iray);
0069 rayinit.ydNpar0 = rayinit.ydNpar0(iray);
0070 rayinit.yP0_2piRp = rayinit.yP0_2piRp(iray);
0071 %
0072 if isfield(rayinit,'yib'),
0073     rayinit.yib = rayinit.yib(iray);
0074 end
0075 %
0076 if isfield(rayinit,'yir'),
0077     rayinit.yir = rayinit.yir(iray);
0078 end
0079 %
0080 if isfield(rayinit,'yim'),
0081     rayinit.yim = rayinit.yim(iray);
0082 end
0083 %
0084 if isfield(rayinit,'yit'),
0085     rayinit.yit = rayinit.yit(iray);
0086 end
0087 %
0088 if isfield(rayinit,'yXfrac'),
0089     rayinit.yXfrac = rayinit.yXfrac(iray);
0090 else
0091     rayinit.yXfrac = NaN;
0092 end
0093 %
0094 tau_lim_init = rayparam.tau_lim;%Initial optical depth limit
0095 kmode = waveparam.kmode;
0096 ydNpar0 = rayinit.ydNpar0;
0097 %
0098 if isfield(rayinit,'launch'),
0099     launch0 = rayinit.launch;
0100 else
0101     launch0 = '';
0102 end
0103 %
0104 if isfield(launch0,'mmode') && launch0.mmode ~= waveparam.mmode,
0105     disp('WARNING : waveparam.mmode enforced with the value from launch0.mmode');
0106     waveparam.mmode = launch0.mmode;
0107 end
0108 %
0109 rays{1}.P0_2piRp = rayinit.yP0_2piRp;%Initial lineic power density in the ray (W/m) [1,1]
0110 %
0111 rays{1}.sx = [];% Ray location (horizontal) along the propagation [1,ns]
0112 rays{1}.sy = [];% Ray location (vertical) along the propagation [1,ns]
0113 rays{1}.sz = [];% Ray location (axial or torodial) along the propagation [1,ns]
0114 %
0115 rays{1}.ss = [];% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0116 %
0117 rays{1}.spsin = [];% Ray location (normalized magnetic flux) along the propagation [1,ns]
0118 rays{1}.stheta = [];% Ray location (poloidal angle) along the propagation [1,ns]
0119 rays{1}.sphi = [];% Ray location (toroidal angle) along the propagation [1,ns]
0120 %
0121 rays{1}.skx = [];% Ray wave vector (major radius) along the propagation [1,ns]
0122 rays{1}.sky = [];% Ray wave vector (vertical) along the propagation [1,ns]
0123 rays{1}.skz = [];% Ray wave vector (axial or toroidal) along the propagation [1,ns]
0124 %
0125 rays{1}.skrho = [];% Ray wave number (radial) along the propagation [1,ns]
0126 rays{1}.sm = [];% Ray wave number (poloidal) along the propagation [1,ns]
0127 rays{1}.sn = [];% Ray wave number (toroidal) along the propagation [1,ns]
0128 %
0129 rays{1}.sNpar = [];% Ray parallel wave number [1,ns]
0130 rays{1}.sNperp = [];% Ray perpendicular wave number [1,ns]
0131 rays{1}.sdNpar = [];% Ray spectral width in Npar [1,ns]
0132 %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. ->
0133 %       square spectrum (centered on sNpar)
0134 rays{1}.sepol_pmz = zeros(3,0);% polarization vector [3,ns] complex
0135 rays{1}.sphi_xyz = zeros(3,0);% normalized power flow [3,ns]
0136 %
0137 rays{1}.spmode = [];%(1) cold mode p, (0) cold mode m (p is the slow LH wave)
0138 %
0139 rays{1}.sr = [];% Ray radial position [1,ns] (m)
0140 rays{1}.srho = [];% Ray normalized radial position [1,ns]
0141 %
0142 rays{1}.salphaphi_lin = [];% Ray absorption coefficient [1,ns]
0143 rays{1}.stau_lin = [];% Ray optical depth [1,ns]
0144 %
0145 rays{1}.seikval = [];% Eikonal approximation validity (OK if less than 1) [1,ns]
0146 rays{1}.sD = [];% Dispersion relation [1,ns]
0147 %
0148 rays{1}.sB = [];% Total magnetic field along the ray [1,ns] (T)
0149 rays{1}.sBRHO = [];% Radial magnetic field along the ray [1,ns] (T)
0150 rays{1}.sBP = [];% Poloidal magnetic field along the ray [1,ns] (T)
0151 rays{1}.sBz = [];% Toroidal or axial magnetic field along the ray [1,ns] (T)
0152 
0153 rays{1}.sTe = [];% Electron temperature along the ray [1,ns] (keV)
0154 rays{1}.sne = [];% Electron density along the ray [1,ns] (m-3)
0155 rays{1}.szTi = zeros(length(equil_fit.zZi),0);% Ion temperatures along the ray [ni,ns] (keV)
0156 rays{1}.szni = zeros(length(equil_fit.zZi),0);% Ion densities along the ray [ni,ns] (m-3)
0157 %
0158 rays{1}.swpe2 = [];%Normalized electron plasma frequency
0159 rays{1}.swce2 = [];%Normalized electron gyro frequency
0160 %
0161 rayinit.tau0_lin = 0;%initial optical depth
0162 rayinit.ss0 = 0;%initial ray length
0163 %
0164 rays{1}.rayinits{1} = rayinit;
0165 %
0166 if isnan(rayinit.yXfrac),% polarization enforced by waveparam.mmode
0167     mmode = sign(waveparam.mmode);%initial polarization
0168 else%
0169     %  polarization enforced by rayinit.yXfrac.
0170     %
0171     mmodethres = 1 - abs(waveparam.mmode);% waveparam.mmode sets the threshold for considering both polarizations.
0172     %
0173     if mmodethres > 1/2,
0174         mmodethres = 1 - mmodethres;% threshold is defined between 0 and 1/2
0175     end
0176     %
0177     if rayinit.yXfrac <= mmodethres,% enforce O mode
0178         mmode = -(1-rayinit.yXfrac);
0179     elseif rayinit.yXfrac >= (1 - mmodethres),% enforce X mode
0180         mmode = rayinit.yXfrac;
0181     else% split power according to rayinit.yXfrac
0182         mmode = [rayinit.yXfrac,1-rayinit.yXfrac];%[Xfrac, Ofrac]
0183     end
0184 end
0185 %
0186 if length(mmode) == 1,% one mode only
0187     rays{1}.rayinits{1}.mmode = sign(mmode);
0188     rays{1}.rayinits{1}.yP0_2piRp = rays{1}.rayinits{1}.yP0_2piRp*abs(mmode);
0189     rays{1}.P0_2piRp = rays{1}.P0_2piRp*abs(mmode);
0190 else% split power according to mmode
0191     %
0192     rays{2} = rays{1};% copy ray parameters
0193     %
0194     rays{1}.rayinits{1}.mmode = 1;% X mode
0195     rays{1}.rayinits{1}.yP0_2piRp = rays{1}.rayinits{1}.yP0_2piRp*mmode(1);
0196     rays{1}.P0_2piRp = rays{1}.P0_2piRp*mmode(1);
0197     %
0198     rays{2}.rayinits{1}.mmode = -1;% O mode
0199     rays{2}.rayinits{1}.yP0_2piRp = rays{2}.rayinits{1}.yP0_2piRp*mmode(2);
0200     rays{2}.P0_2piRp = rays{2}.P0_2piRp*mmode(2);
0201     %
0202 end  
0203 %
0204 % --------------
0205 %
0206 % rayparam.reflection :
0207 %     * 0 : no reflection
0208 %    * 1 : reflections in C3PO on LCFS (LH waves...)
0209 %    * -n : up to n reflections on vessel internal wall
0210 %   * (+1i) : calc rayinit for last reflection, no new C3PO calculation (ex: stray emission)
0211 %
0212 lastref = imag(rayparam.reflection);% calc rayinit for last reflection, no new C3PO calculation (ex: stray emission)
0213 rayparam.reflection = real(rayparam.reflection);
0214 %
0215 nref = max([0,-rayparam.reflection]);
0216 rayparam.reflection = max([rayparam.reflection,0]);% take only positive or null in C3PO
0217 %
0218 if nref > 0,
0219     C3POdisplay.ray = 0;
0220 end
0221 %
0222 for iref = 1:nref + 1,% vessel reflection index
0223     %
0224     nrays = length(rays);
0225     oldrays = zeros(1,nrays);
0226     newrays = cell(1,0);
0227     %
0228     for irayref = 1:nrays,% each ray generated by vessel reflection must be considered separately
0229         %
0230         if length(rays{irayref}.rayinits) < iref,% this ray is absorbed or not reflected
0231             continue
0232         end
0233         %
0234         rayinit = rays{irayref}.rayinits{end};
0235         %
0236         y0 = init_RT_jd(equil_fit,...
0237             rayinit.yrho0,rayinit.ytheta0,rayinit.yz0,omega_rf,...
0238             rayinit.yNpar0,rayinit.ym0,rayinit.ykz0,rayinit.mmode,kmode,fluct_fit).';
0239         %
0240         if C3POdisplay.ray,
0241             info_dke_yp(2,['Initial conditions of ray #',int2str(iray),' are calculated']);
0242         end
0243         %
0244         if isreal(y0(4)) && abs(y0(4)) > 0,%ray enters the plasma
0245             %
0246             rayparam.tau_lim = tau_lim_init - rayinit.tau0_lin;% the absorption criterion accounts for previous steps
0247             %
0248             [tout,yout,yproc] = raytracing_yp(rayparam,equil_fit,y0,f0struct,fluct_fit,flag);
0249             %
0250             if strcmp(computer,'MACI64'),
0251                 clear mex;%contribute to stabilize the ray-tracing. At least for mac...experimental at this stage
0252             end
0253             %
0254             flag_rhomin = false;
0255             %
0256             while ~isempty(yout) && yout(end,1) <= rayparam.rhomin,
0257                 %
0258                 flag_rhomin = true;
0259                 %
0260                 if ~isinf(equil_fit.Rp),%torus
0261                     %
0262                     X0 = (yproc(end-2,14)*equil_fit.ap + equil_fit.Rp).*cos(yout(end-2,3)/equil_fit.Rp);
0263                     Y0 = - (yproc(end-2,14)*equil_fit.ap  + equil_fit.Rp).*sin(yout(end-2,3)/equil_fit.Rp);%Note : phi is clockwise from top
0264                     Z0 = yproc(end-2,15)*equil_fit.ap + equil_fit.Zp;
0265                     %
0266                     X1 = (yproc(end-1,14)*equil_fit.ap + equil_fit.Rp).*cos(yout(end-1,3)/equil_fit.Rp);%input point
0267                     Y1 = - (yproc(end-1,14)*equil_fit.ap  + equil_fit.Rp).*sin(yout(end-1,3)/equil_fit.Rp);%input point
0268                     Z1 = yproc(end-1,15)*equil_fit.ap + equil_fit.Zp;%input point
0269                     %
0270                     d = sqrt((X1 - X0)^2 + (Y1 - Y0)^2 + (Z1 - Z0)^2);
0271                     %
0272                     Vx = (X1 - X0)/d;%local direction of the ray
0273                     Vy = (Y1 - Y0)/d;%local direction of the ray
0274                     Vz = (Z1 - Z0)/d;%local direction of the ray
0275                     V = [Vx,Vy,Vz];
0276                     %
0277                     R_ref = sqrt(X1.^2 + Y1.^2);
0278                     Z_ref = Z1;
0279                     phi_ref = -atan(Y1./X1) + pi*(X1 < 0) - 2*pi*(X1 < 0 && Y1 > 0);% -pi < phi < pi
0280                     %
0281                     NR_ref = Vx*cos(phi_ref) - Vy*sin(phi_ref);
0282                     Nphi_ref = -Vx*sin(phi_ref) - Vy*cos(phi_ref);
0283                     NZ_ref = Vz;
0284                     %
0285                     launch1 = launch0;
0286                     launch1.launch0 = launch0;
0287                     %
0288                     launch1.yR_L = R_ref;
0289                     launch1.yZ_L = Z_ref;
0290                     launch1.yphi_L = phi_ref;
0291                     launch1.yalpha_L = - atan(Nphi_ref/NR_ref) + pi*(NR_ref < 0) - 2*pi*(NR_ref < 0 && Nphi_ref > 0);% -pi < alpha < pi
0292                     launch1.ybeta_L = acos(NZ_ref);% 0 < beta < pi
0293                     %
0294                     [rayinit_guess,yout1] = main_rayinit_launch_jd(equil,launch1,rayparam.rhomin);
0295                     %
0296                     rayinit_guess.yNpar0 = (yproc(end-1,2) - yproc(end-2,2))/(yout(end-1,end) - yout(end-2,end))*max(yout1(:,end))+ yproc(end-1,2);%linear npar extrapolation for a continuity of Npar.
0297                     %
0298                     yout1(:,6) = yout(1,6);%axisymmetry. kz is constant
0299                     %
0300                     % start m optimization
0301                     %
0302                     rayparam_guess = rayparam;
0303                     %
0304                     rayparam_guess.tfinal = 300;%not too smal to avoid memory problem in C3PO
0305                     %
0306                     rayinit_guess.tau0_lin = rayinit.tau0_lin;
0307                     rayinit_guess.ss0 = rayinit.ss0;
0308                     rayinit_guess.mmode = rayinit.mmode + j;%j is for outward propagation from the considered plasma domain
0309                     %
0310                     m_guess = (yout(end,5) + yout1(end,5))/2;%mean value between the input and output
0311                     %
0312                     %delta_guess = mguess_yp(equil_fit,f0struct,fluct_fit,flag,kmode,rayparam_guess,rayinit_guess,V,m_guess)
0313                     %
0314                     [m_opt,delta_opt] = fminsearch(@(m_guess) mguess_yp(equil_fit,f0struct,fluct_fit,flag,kmode,rayparam_guess,rayinit_guess,V,m_guess),m_guess)
0315                     %
0316                     % propagation inside the torus
0317                     %
0318                     yout1(:,5) = yout1(:,end)*(m_opt - yout(end-1,5))/max(yout1(:,end)) + yout(end-1,5);%progressive evolution of m inside the torus
0319                     Npar_in = yout1(:,end)*(rayinit_guess.yNpar0 - yproc(end-1,2))/max(yout1(:,end)) + yproc(end-1,2);%progressive evolution of Npar inside the torus
0320                     %
0321                     yout1a = [];
0322                     %
0323                     for ii = 1:length(yout1)
0324                         yout1a = [yout1a;[init_RT_jd(equil_fit,...
0325                         yout1(ii,1),yout1(ii,2),yout1(ii,3),omega_rf,...
0326                         Npar_in(ii),yout1(ii,5),NaN,rayinit_guess.mmode,kmode,fluct_fit).',NaN]];
0327                     end
0328                     %
0329                     tout = [tout(1:end-1,:);zeros(length(yout1a),1)];
0330                     yout = [yout(1:end-1,:);yout1a];
0331                     %
0332                 else% cylindrical case
0333                     %
0334                     disp('WARNING: cylindrical case not yet implemented');
0335                     keyboard
0336                 end
0337                 %
0338                 % propagation after the torus
0339                 %
0340                 rayinit1.tau0_lin = rayinit.tau0_lin;
0341                 rayinit1.ss0 = rayinit.ss0;
0342                 rayinit1.mmode = rayinit.mmode + j;%j is for outward propagation from the considered plasma domain
0343                 %
0344                 y1 = init_RT_jd(equil_fit,...
0345                     yout(end,1),yout(end,2),yout(end,3),omega_rf,...
0346                     rayinit_guess.yNpar0,yout(end,5),NaN,rayinit_guess.mmode,kmode,fluct_fit).';
0347                 %
0348                 [tout1, yout1, yproc1] = raytracing_yp(rayparam,equil_fit,y1,f0struct,fluct_fit,flag);
0349                 %
0350                 tout = [tout;tout1];
0351                 yout = [yout;yout1];
0352                 %
0353             end
0354             %
0355             if flag_rhomin,
0356                 yproc = [];%ray properties are not calculated in the C program is rhomin is catched
0357             end
0358             %
0359             if C3POdisplay.ray,
0360                 info_dke_yp(2,['Trajectory of ray #',int2str(iray),' is calculated']);
0361             end
0362             %
0363             filename_in = [id_wave,'_ray_',int2str(iray)];
0364             %
0365             [iabs_out,srho,stheta,sx,sy,sz,sNpar,sNperp,sepol_pmz,sphi_xyz,spmode,sD,ss,salphaphi_lin,sB,sBRHO,sBP,sBz,sTe,sne,szTi,szni,sr,salpha,sgradrhon,stau_lin,seikval,skrho,sm,skz,sn,swpe2,swce2] = raypostprocessing_yp(tout,yout,yproc,rayparam,equil_fit,fluct_fit,C3POdisplay.text,C3POdisplay.ray,C3POdisplay.p_opt,filename_in);
0366             %
0367             if C3POdisplay.ray == 1,
0368                 pause(5)
0369             end
0370             %
0371             % building ray for FP calculations
0372             %
0373             rays{irayref}.sx = [rays{irayref}.sx,sx.'];% Ray location (horizontal) along the propagation [1,ns]
0374             rays{irayref}.sy = [rays{irayref}.sy,sy.'];% Ray location (vertical) along the propagation [1,ns]
0375             rays{irayref}.sz = [rays{irayref}.sz,sz.'];% Ray location (axial) along the propagation [1,ns]
0376             %
0377             rays{irayref}.ss = [rays{irayref}.ss,rayinit.ss0 + ss.'];% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0378             rays{irayref}.stheta = [rays{irayref}.stheta,stheta.'];% Ray location (poloidal angle) along the propagation [1,ns]
0379             %
0380             if ~isinf(equil_fit.Rp),
0381                 sphi = sz/equil_fit.Rp;
0382             else
0383                 sphi = NaN(size(sz));
0384             end
0385             %
0386             rays{irayref}.sphi = [rays{irayref}.sphi,sphi.'];
0387             %
0388             rays{irayref}.skrho = [rays{irayref}.skrho,skrho.'];% Ray radial wave vector [1,ns]
0389             rays{irayref}.sm = [rays{irayref}.sm,sm.'];% Ray poloidal wave number [1,ns]
0390             rays{irayref}.skz = [rays{irayref}.skz,skz.'];% Ray axial wave vector [1,ns]
0391             rays{irayref}.sn = [rays{irayref}.sn,sn.'];% Ray toroidal wave number [1,ns]
0392             %
0393             rays{irayref}.sNpar = [rays{irayref}.sNpar,sNpar.'];% Ray parallel wave number [1,ns]
0394             rays{irayref}.sNperp = [rays{irayref}.sNperp,sNperp.'];% Ray perpendicular wave number [1,ns]
0395             rays{irayref}.sdNpar = [rays{irayref}.sdNpar,ydNpar0*ones(1,length(srho))];% Ray spectral width in Npar [1,ns]
0396             %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. ->
0397             %       square spectrum (centered on sNpar)
0398             rays{irayref}.sepol_pmz = [rays{irayref}.sepol_pmz,sepol_pmz.'];% polarization vector [3,ns] complex
0399             rays{irayref}.sphi_xyz = [rays{irayref}.sphi_xyz,sphi_xyz.'];% normalized power flow [3,ns]
0400             %
0401             rays{irayref}.spmode = [rays{irayref}.spmode,spmode.'];%(1) cold mode p, (0) cold mode m (p is the slow LH wave)
0402             %
0403             rays{irayref}.sr = [rays{irayref}.sr,sr.'];% Ray radial position [1,ns] (m)
0404             rays{irayref}.srho = [rays{irayref}.srho,srho.'];% Ray normalized radial position [1,ns]
0405             %
0406             rays{irayref}.salphaphi_lin = [rays{irayref}.salphaphi_lin,salphaphi_lin.'];% Ray absorption coefficient [1,ns]
0407             rays{irayref}.stau_lin = [rays{irayref}.stau_lin,rayinit.tau0_lin + stau_lin.'];% Ray optical depth [1,ns]
0408             %
0409             rays{irayref}.seikval = [rays{irayref}.seikval,seikval.'];% Eikonal approximation validity (OK if less than 1) [1,ns]
0410             rays{irayref}.sD = [rays{irayref}.sD,sD.'];% Dispersion relation [1,ns]
0411             %
0412             rays{irayref}.sB = [rays{irayref}.sB,sB.'];% Total magnetic field along the ray [1,ns] (T)
0413             rays{irayref}.sBRHO = [rays{irayref}.sBRHO,sBRHO.'];% Radial magnetic field along the ray [1,ns] (T)
0414             rays{irayref}.sBP = [rays{irayref}.sBP,sBP.'];% Poloidal magnetic field along the ray [1,ns] (T)
0415             rays{irayref}.sBz = [rays{irayref}.sBz,sBz.'];% Toroidal magnetic field along the ray [1,ns] (T)
0416             
0417             rays{irayref}.sTe = [rays{irayref}.sTe,sTe.'];% Electron temperature along the ray [1,ns] (keV)
0418             rays{irayref}.sne = [rays{irayref}.sne,sne.'];% Electron density along the ray [1,ns] (m-3)
0419             rays{irayref}.szTi = [rays{irayref}.szTi,szTi.'];% Ion temperatures along the ray [ni,ns] (keV)
0420             rays{irayref}.szni = [rays{irayref}.szni,szni.'];% Ion densities along the ray [ni,ns] (m-3)
0421             %
0422             rays{irayref}.swpe2 = [rays{irayref}.swpe2,swpe2];%Normalized electron plasma frequency
0423             rays{irayref}.swce2 = [rays{irayref}.swce2,swce2];%Normalized electron gyro frequency
0424             %
0425             if (iref <= nref || lastref) && yout(end,1) > 1 - 2*rayparam.dS && ~isinf(equil.Rp),% ray reaches edge before full absorption and a reflection is enforced
0426                 %
0427                 sR = sx + equil.Rp;
0428                 sZ = sy + equil.Zp;
0429                 %
0430                 R_end = sR(end); 
0431                 phi_end = sphi(end);
0432                 Z_end = sZ(end);        
0433                 %
0434                 X_end = R_end.*cos(phi_end);
0435                 Y_end = -R_end.*sin(phi_end);%note phi is of opposite sign to usual convention (it is clockwise from top)
0436                 %
0437                 krho_end = yout(end,4);
0438                 m_end = yout(end,5);
0439                 n_end = yout(end,6);
0440                 %
0441                 gradrhon_end = sgradrhon(end);
0442                 alpha_end = salpha(end);
0443                 theta_end = stheta(end);
0444                 r_end = sr(end);
0445                 %
0446                 kR_end = krho_end*gradrhon_end*cos(theta_end - alpha_end)/equil_fit.ap - m_end*sin(theta_end)/r_end;
0447                 kZ_end = krho_end*gradrhon_end*sin(theta_end - alpha_end)/equil_fit.ap + m_end*cos(theta_end)/r_end;
0448                 kphi_end = n_end/R_end;
0449                 %
0450                 kX_end = kR_end*cos(phi_end) - kphi_end*sin(phi_end);
0451                 kY_end = -kR_end*sin(phi_end) - kphi_end*cos(phi_end);
0452                 %
0453                 P_end = [X_end;Y_end;Z_end];
0454                 K_end = [kX_end;kY_end;kZ_end];
0455                 %
0456                 % calculate the polarization in the zero-density limit so that K is perpendicular to E
0457                 %
0458                 B_end = sB(end);
0459                 BPHI_end = sBz(end);
0460                 Bx_end = sBP(end)*sin(alpha_end - theta_end);%note that sBP is signed
0461                 By_end = sBP(end)*cos(alpha_end - theta_end);
0462                 %
0463                 BX_end = Bx_end*cos(phi_end) - BPHI_end*sin(phi_end);
0464                 BY_end = -Bx_end*sin(phi_end) - BPHI_end*cos(phi_end);
0465                 BZ_end = By_end;
0466                 %
0467                 B_end_XYZ = [BX_end;BY_end;BZ_end];
0468                 %
0469                 N_end = K_end/norm(K_end);%index of refraction is 1 in vacuum
0470                 %
0471                 z_XYZ = B_end_XYZ/B_end;
0472                 y_XYZ = cross(B_end_XYZ,N_end)/norm(cross(B_end_XYZ,N_end));
0473                 x_XYZ = cross(y_XYZ,z_XYZ);
0474                 %
0475                 M_xyzXYZ = [x_XYZ,y_XYZ,z_XYZ];
0476                 %
0477                 Npar_end = N_end.'*z_XYZ;%value in vacuum
0478                 %
0479                 [Nperpp,Nperpm,Kxyz,ep_xyz,em_xyz] = colddisp_dke_jd(Npar_end,omega_rf,NaN,0,1,1,B_end);
0480                 %
0481                 ep_XYZ = M_xyzXYZ*ep_xyz;
0482                 em_XYZ = M_xyzXYZ*em_xyz;
0483                 %
0484                 if spmode(end) == 1,
0485                     E_end = ep_XYZ;
0486                 else
0487                     E_end = em_XYZ;
0488                 end
0489                 %
0490                 %E_end = sepol_pmz(end,:).';
0491                 %
0492                 nrefvac = 10;% number of vacuum reflections allowed
0493                 launch_back = '';
0494                 %
0495                 for i=1:nrefvac% to allow reflections without reentering the plasma
0496                     %
0497                     %try
0498                         save rayend.mat P_end N_end E_end
0499                         %keyboard
0500                         [P_ref,N_ref,E_ref] = reflection_jd(equil,P_end,N_end,E_end);
0501                     %catch
0502                         %keyboard
0503                     %end
0504                     %
0505                     R_ref = sqrt(P_ref(1).^2 + P_ref(2).^2);
0506                     Z_ref = P_ref(3);
0507                     phi_ref = -atan(P_ref(2)./P_ref(1)) + pi*(P_ref(1) < 0) - 2*pi*(P_ref(1) < 0 && P_ref(2) > 0);% -pi < phi < pi
0508                     %
0509                     NR_ref = N_ref(1)*cos(phi_ref) - N_ref(2)*sin(phi_ref);
0510                     Nphi_ref = -N_ref(1)*sin(phi_ref) - N_ref(2)*cos(phi_ref);
0511                     NZ_ref = N_ref(3);
0512                     %
0513                     launch = launch0;
0514                     launch.launch_back = launch_back;
0515                     %
0516                     launch.yR_L = R_ref;
0517                     launch.yZ_L = Z_ref;
0518                     launch.yphi_L = phi_ref;
0519                     launch.yalpha_L = - atan(Nphi_ref/NR_ref) + pi*(NR_ref < 0) - 2*pi*(NR_ref < 0 && Nphi_ref > 0);% -pi < alpha < pi
0520                     launch.ybeta_L = acos(NZ_ref);% 0 < beta < pi
0521                     %
0522                     launch.N_XYZ = N_ref;
0523                     launch.epol_XYZ = E_ref;
0524                     %
0525                     rayinit = main_rayinit_launch_jd(equil,launch);
0526                     %
0527                     if ~isnan(rayinit.yrho0),% the ray reenters the plasma
0528                         break
0529                     else
0530                         launch_back = launch;
0531                         P_end = P_ref;
0532                         N_end = N_ref;
0533                         E_end = E_ref;
0534                     end
0535                 end
0536                 %
0537                 rayinit.tau0_lin = rays{irayref}.stau_lin(end);
0538                 rayinit.ss0 = rays{irayref}.ss(end);
0539                 %
0540                 if iref == nref + 1,% end of calculation, just add last point of reflection,
0541                     %
0542                     rays{irayref}.rayinits = [rays{irayref}.rayinits,rayinit];
0543                     %
0544                 elseif ~isempty(rayinit.yrho0) && ~isnan(rayinit.yrho0),% The reflected ray reaches the plasma, new rays are formed
0545                     %
0546                     ray = rays{irayref};
0547                     %
0548                     rayX = ray;
0549                     rayO = ray;
0550                     %
0551                     rayX.P0_2piRp = ray.P0_2piRp*rayinit.am2;
0552                     %
0553                     rayO.P0_2piRp = ray.P0_2piRp*rayinit.ap2;
0554                     %
0555                     rayinit.mmode = 1;
0556                     rayX.rayinits = [rayX.rayinits,rayinit];
0557                     %
0558                     rayinit.mmode = -1;
0559                     rayO.rayinits = [rayO.rayinits,rayinit];
0560                     %
0561                     newrays = [newrays,rayX,rayO];% two new rays obtained from duplicating original ray
0562                     oldrays(irayref) = 1;% old duplicated ray to be removyed
0563                     %
0564                 end
0565             end
0566         else
0567             disp('Pb krho at the edge. No ray tracing');
0568         end
0569     end
0570     %
0571     rays = [rays(oldrays == 0),newrays];
0572     %
0573 end
0574 %
0575 if nargout > 1,
0576    %
0577    rhomin = 1;
0578    %
0579    for iray = 1:length(rays),
0580        rhomin = min([rhomin,rays{iray}.srho]);
0581    end
0582 end
0583 end
0584 %
0585 function [delta] = calc_rhomin_jd(ym0,iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,waveparam,C3POdisplay,f0struct,fluct_fit)
0586     %
0587     ym0_init = rayinit.ym0;
0588     %
0589     rayinit.ym0 = ym0;
0590     %
0591     [~,rhomin] = loop_main_C3PO_jd(iray,equil,equil_fit,id_wave,omega_rf,rayinit,rayparam,flag,rmfield(waveparam,'rhomin'),C3POdisplay,f0struct,fluct_fit);
0592     %
0593     delta = abs(waveparam.rhomin - rhomin);
0594     %
0595     disp(' ');
0596     disp([' ym0 : ',num2str(ym0_init),' -> ',num2str(ym0),' ; rhomin : ',num2str(rhomin),' ; delta : ',num2str(delta)]);
0597     %
0598 end
0599 
0600 function delta = mguess_yp(equil_fit,f0struct,fluct_fit,flag,kmode,rayparam_guess,rayinit_guess,V,m_guess)
0601     %
0602     y0_guess = init_RT_jd(equil_fit,...
0603     rayinit_guess.yrho0,rayinit_guess.ytheta0,rayinit_guess.yz0,rayinit_guess.omega_rf,...
0604     rayinit_guess.yNpar0,m_guess,rayinit_guess.ykz0,rayinit_guess.mmode,kmode,fluct_fit).';
0605     %
0606     [tout_guess, yout_guess, yproc_guess] = raytracing_yp(rayparam_guess,equil_fit,y0_guess,f0struct,fluct_fit,flag);
0607     %
0608     X2 = (yproc_guess(1,14)*equil_fit.ap + equil_fit.Rp).*cos(yout_guess(1,3)/equil_fit.Rp);
0609     Y2 = - (yproc_guess(1,14)*equil_fit.ap  + equil_fit.Rp).*sin(yout_guess(1,3)/equil_fit.Rp);%Note : phi is clockwise from top
0610     Z2 = yproc_guess(1,15)*equil_fit.ap + equil_fit.Zp;
0611     %
0612     X3 = (yproc_guess(2,14)*equil_fit.ap + equil_fit.Rp).*cos(yout_guess(2,3)/equil_fit.Rp);%input point
0613     Y3 = - (yproc_guess(2,14)*equil_fit.ap  + equil_fit.Rp).*sin(yout_guess(2,3)/equil_fit.Rp);%input point
0614     Z3 = yproc_guess(2,15)*equil_fit.ap + equil_fit.Zp;%input point
0615     %
0616     d_guess = sqrt((X3 - X2)^2 + (Y3 - Y2)^2 + (Z3 - Z2)^2);
0617     %
0618     Vx_guess = (X3 - X2)/d_guess;%local direction of the ray
0619     Vy_guess = (Y3 - Y2)/d_guess;%local direction of the ray
0620     Vz_guess = (Z3 - Z2)/d_guess;%local direction of the ray
0621     V_guess = [Vx_guess,Vy_guess,Vz_guess];
0622     %
0623     delta = 1 - V(:)'*V_guess(:);
0624 end
0625