main_rayinit_launch_jd

 This function creates the rayinit structure based on ray launching parameters. Straight ray
 propagation is assumed up to the last closed flux-surface for the EC wave

 INPUTS:
   Launch structure fields:

       - type = Type of wave (LH or EC;
       - omega_rf = [4.6]*2*pi*1e9;%wave angular frequency (rad/s)

   Specific to LH wave

       - bPlhtot_2piRp: initial lineic power density in the ray (W/m) [1,n_lobe]
       - bNpar0: Initial Npar value [1,n_lobe]
       - bdNpar0: Initial dNpar value [1,n_lobe]
       - ry0: vertical launching positions [m_pol,1]

       - m0: initial poloidal mode number
       - kz0: initial axial wave vector
       - z0: initial axial position
       - i_ref: Initial radial position index
       - LFS: (-1) high field side launching, (1) low field side launching

   Specific to EC wave

       - yx_L : major radius position of launching (m)
       - yy_L : vertical position of launching (m)
       - yz_L : axial position of launching (m), (R,Z,phi) or (x,y,z) right handed
       - yalpha_L : horizontal angle with respect to R (rad)
       - ybeta_L : vertical angle with respect to Z (rad) 
       - yP_L_2piRp : lineic power density launched in the ray (W/m)
       - dNpar0 :  initial power spectrum width. May be NaN
       - ns : (1000)
       - method1: 'pchip' for 1D interpolation
        - method2: 'cubic' for 2D interpolation

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

0001 function [rayinit,yout] = main_rayinit_launch_jd(equil,launch,rhomin)
0002 %
0003 % This function creates the rayinit structure based on ray launching parameters. Straight ray
0004 % propagation is assumed up to the last closed flux-surface for the EC wave
0005 %
0006 % INPUTS:
0007 %   Launch structure fields:
0008 %
0009 %       - type = Type of wave (LH or EC;
0010 %       - omega_rf = [4.6]*2*pi*1e9;%wave angular frequency (rad/s)
0011 %
0012 %   Specific to LH wave
0013 %
0014 %       - bPlhtot_2piRp: initial lineic power density in the ray (W/m) [1,n_lobe]
0015 %       - bNpar0: Initial Npar value [1,n_lobe]
0016 %       - bdNpar0: Initial dNpar value [1,n_lobe]
0017 %       - ry0: vertical launching positions [m_pol,1]
0018 %
0019 %       - m0: initial poloidal mode number
0020 %       - kz0: initial axial wave vector
0021 %       - z0: initial axial position
0022 %       - i_ref: Initial radial position index
0023 %       - LFS: (-1) high field side launching, (1) low field side launching
0024 %
0025 %   Specific to EC wave
0026 %
0027 %       - yx_L : major radius position of launching (m)
0028 %       - yy_L : vertical position of launching (m)
0029 %       - yz_L : axial position of launching (m), (R,Z,phi) or (x,y,z) right handed
0030 %       - yalpha_L : horizontal angle with respect to R (rad)
0031 %       - ybeta_L : vertical angle with respect to Z (rad)
0032 %       - yP_L_2piRp : lineic power density launched in the ray (W/m)
0033 %       - dNpar0 :  initial power spectrum width. May be NaN
0034 %       - ns : (1000)
0035 %       - method1: 'pchip' for 1D interpolation
0036 %        - method2: 'cubic' for 2D interpolation
0037 %
0038 % by J. Decker (CEA-IRFM) <joan.decker@cea.fr> and Y. Peysson (IRFM/DSM/CEA) <yves.peysson@cea.fr>
0039 %
0040 if nargin < 2,
0041     error('Not enough input arguments in main_rayinit_launch_jd');
0042 end
0043 %
0044 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0045 %
0046 if ~isfield(launch,'m0'),
0047     launch.m0 = NaN;
0048 end
0049 %
0050 if ~isfield(launch,'n0'),
0051     launch.n0 = NaN;
0052 end
0053 %
0054 if ~isfield(launch,'kz0'),
0055     launch.kz0 = NaN;
0056 end
0057 %
0058 if strcmp(launch.type,'EC') && ~isfield(launch,'dNpar0'),
0059     launch.dNpar0 = NaN;
0060 end
0061 if ~isfield(launch,'ns'),
0062     launch.ns = 1000;
0063 end
0064 if ~isfield(launch,'method1'),
0065     launch.method1 = 'pchip';
0066 end
0067 %
0068 if ~isfield(launch,'method2'),
0069     launch.method2 = 'spline';
0070 end
0071 %
0072 if strcmp(launch.type,'EC'),
0073     %
0074     dNpar0 = launch.dNpar0;
0075     ns = launch.ns;
0076     method1 = launch.method1;
0077     method2 = launch.method2;
0078     %
0079     if isfield(launch,'yP_L') && ~isinf(equil.Rp),
0080         yP0_2piRp = launch.yP_L/(2*pi*equil.Rp);%initial power in the ray (W) [1,n_lobe]
0081     elseif isfield(launch,'yP_L_2piRp')
0082         yP0_2piRp = launch.yP_L_2piRp;
0083     else
0084         error('ERROR in main_rayinit_launch_jd.m -> inconsistent launch input: use yP_L if Rp is finite only, or yP_L_2piRp otherwise.');
0085     end
0086     %
0087     ny = length(yP0_2piRp);
0088     %
0089     ptr = sqrt(equil.ptx.^2 + equil.pty.^2);
0090     smax = 2*max(max(ptr));
0091     s = linspace(0,smax,ns);
0092     %
0093     % Straight ray propagation to reach first closed flux surface
0094     %
0095     if isfield(launch,'yR_L') && ~isinf(equil.Rp),
0096         yx_L = launch.yR_L - equil.Rp;
0097     elseif isfield(launch,'yx_L')
0098         yx_L = launch.yx_L;
0099     else
0100         error('ERROR in main_rayinit_launch_jd.m -> inconsistent launch input: use yR_L if Rp is finite only, or yx_L otherwise.');
0101     end
0102     %
0103     if isfield(launch,'yZ_L')
0104         yy_L = launch.yZ_L - equil.Zp;
0105     elseif isfield(launch,'yy_L')
0106         yy_L = launch.yy_L;
0107     end
0108     %
0109     if isfield(launch,'yphi_L') && ~isinf(equil.Rp),
0110         yz_L = launch.yphi_L*equil.Rp;
0111     elseif isfield(launch,'yz_L')
0112         yz_L = launch.yz_L;
0113     else
0114         error('ERROR in main_rayinit_launch_jd.m -> inconsistent launch input: use yphi_L if Rp is finite only, or yz_L otherwise.');
0115     end
0116     %
0117     [sysh,sysy,syx,syy,syz,syalpha,sybeta,sygx,sygy,sygz] = raypath_prop_jd(equil.Rp,yx_L,yy_L,-yz_L,launch.yalpha_L,launch.ybeta_L,s);    
0118     %
0119     % The definition of z in raypath_prop_jd is opposite from the definition in DKE
0120     %
0121     syz = - syz;
0122     sygz = - sygz;
0123     %
0124     syr = sqrt(syx.^2 + syy.^2);
0125     syx(syx == 0) = eps;
0126     sytheta = atan(syy./syx) + pi*(syx < 0) + 2*pi*(syx > 0 & syy < 0);
0127     %
0128     % In vacuum N = g;
0129     %
0130     syNx = sygx;
0131     syNy = sygy;
0132     syNz = sygz;
0133     %
0134     psi_apRp = equil.psi_apRp;
0135     psin = psi_apRp/psi_apRp(end);
0136     %
0137     rayinit.omega_rf = launch.omega_rf;
0138     %
0139     for iy = 1:ny,    
0140         %
0141         if ~isnan(yx_L(iy)),
0142             %
0143             %
0144             % interpolation for r at ray poloidal locations
0145             %
0146             pstr = interp1(equil.theta,ptr.',sytheta(:,iy),method1).';%r(psi,st)
0147             %
0148             % interpolation for psi at ray position (along cst theta lines)
0149             %
0150             %spsin = interp1(pstr,ppsin,sr,method1);%psi(s)
0151             %
0152             spsin = zeros(1,ns);
0153             %
0154             for is = 1:ns 
0155                 spsin(is) = interp1(pstr(:,is),psin,syr(is,iy),method1,NaN);%psi(s)
0156             end
0157             %
0158             if nargin == 2, %ray that enters the plasma
0159                 is0 = find(~isnan(spsin),1,'first');
0160                 is0 = find(spsin <= 1,1,'first');
0161             else %ray that leaves the small torus around the magnetic axis
0162                 is0 = find(spsin >= rho2psi_jd(equil,rhomin)/psi_apRp(end),1,'first');
0163             end
0164             %
0165             psin0 = spsin(is0);
0166             %
0167             s0 = sqrt(sysh(is0,iy).^2 + sysy(is0,iy).^2);%distance along the ray from mirror to LCFS
0168             %
0169             x0 = syx(is0,iy);
0170             y0 = syy(is0,iy);
0171             r0 = syr(is0,iy);
0172             theta0 = sytheta(is0,iy);
0173             z0 = syz(is0,iy);
0174             Nx0 = syNx(is0,iy);
0175             Ny0 = syNy(is0,iy);
0176             Nz0 = syNz(is0,iy);
0177             %
0178             rho0 = psi2rho_jd(equil,psin0);
0179             %
0180         else
0181             %
0182             s0 = 0;%distance along the ray from mirror to LCFS
0183             %
0184             rho0 = 1;
0185             psin0 = 1;
0186             %
0187             tx = equil.ptx(end,:);
0188             ty = equil.pty(end,:);
0189             %
0190             mask = [find(ty(1:end - 1) < 0 & diff(ty) > 0),1 + [0,find(ty(1:end - 1) >= 0 & diff(ty) > 0)]]; 
0191             %
0192             y0 = yy_L(iy);
0193             x0 = interp1(ty(mask),tx(mask),y0,method1,NaN);
0194             %
0195             if isnan(y0),
0196                 error('The vertical launching position in beyond the LCFS');
0197             end
0198             %
0199             r0 = sqrt(x0.^2 + y0.^2);
0200             theta0 = atan(y0./x0) + pi*(x0 < 0) + 2*pi*(x0 > 0 & y0 < 0);
0201             z0 = - yz_L(iy);
0202             %
0203             Nx0 = sin(launch.ybeta_L(iy)).*cos(launch.yalpha_L(iy));
0204             Ny0 = cos(launch.ybeta_L(iy));
0205             Nz0 = sin(launch.ybeta_L(iy)).*sin(launch.yalpha_L(iy));
0206             %
0207         end
0208         %
0209         if nargout == 2 && ny == 1,%for the rhomin > 0 case (EC wave only)
0210             %
0211             yout(1:is0,1) = psi2rho_jd(equil,spsin(1:is0));%rho
0212             yout(1:is0,2) = sytheta(1:is0);%theta
0213             yout(1:is0,3) = syz(1:is0);%z
0214             %
0215             yout(1:is0,4) = NaN;%krho
0216             yout(1:is0,5) = NaN;%m
0217             yout(1:is0,6) = NaN;%kz
0218             %
0219             yout(1:is0,7) = launch.omega_rf;
0220             yout(1:is0,8) = sqrt(sysh(1:is0).^2 + sysy(1:is0).^2);%distance along the ray
0221             %
0222         else
0223             yout = [];
0224         end
0225         %
0226         if ~isempty(psin0),
0227             %
0228             Bx0 = interp2(psin,equil.theta,equil.ptBx.',psin0,theta0,method2);
0229             By0 = interp2(psin,equil.theta,equil.ptBy.',psin0,theta0,method2);
0230             Bz0 = interp2(psin,equil.theta,equil.ptBPHI.',psin0,theta0,method2);    
0231             %
0232             B0 = sqrt(Bx0.^2 + By0.^2 + Bz0.^2);
0233             Npar0 = (Nx0*Bx0 + Ny0*By0 + Nz0*Bz0)/B0;
0234             %
0235             Bp0 = sqrt(Bx0.^2 + By0.^2);
0236             cosalpha0 = (By0*x0 - Bx0*y0)/r0/Bp0;
0237             Ns0 = (Nx0*Bx0 + Ny0*By0)/Bp0;
0238             m0 = r0*Ns0*rayinit.omega_rf/clum/cosalpha0;
0239             %
0240             rayinit.yrho0(iy) = rho0;%Initial radial position at launch
0241             rayinit.ytheta0(iy) = theta0;%Initial poloidal position at launch
0242             rayinit.yz0(iy) = z0;%Initial axial position at launch
0243             rayinit.ym0(iy) = m0;%Initial poloidal mode number
0244             rayinit.ykz0(iy) = launch.kz0;%Initial toroidal mode number
0245             rayinit.yNpar0(iy) = Npar0;%Initial index of refraction
0246             rayinit.ydNpar0(iy) = dNpar0;%initial Ray spectral width
0247             rayinit.yP0_2piRp(iy) = yP0_2piRp(iy);%Lineic initial power density initial power in the ray (W/m)
0248             %
0249             if nargout == 2 && ny == 1,%for the rhomin > 0 case (EC wave only)
0250                 yout(is0,5) = m0;%m0 guess
0251             end
0252             %
0253             rayinit.yib(iy) = iy;% beam index for tracking (fluctuations)
0254             %
0255             if ny == 1 && isfield(launch,'w0'),%beam waist
0256                 rayinit.w0 = launch.w0;
0257                 if isfield(launch,'z_L'),%for divergence calculation
0258                     rayinit.z0 = launch.z_L - s0;%distance from LCFS to waist
0259                 else
0260                     rayinit.z0 = 0;
0261                 end
0262             else% the Gaussian beam representation is compatible with the one ray description only
0263                 rayinit.w0 = NaN;
0264                 rayinit.z0 = NaN;
0265             end
0266             %
0267             if isfield(launch,'yXfrac'),
0268                 rayinit.yXfrac(iy) = launch.yXfrac(iy);
0269             else
0270                 rayinit.yXfrac(iy) = NaN;
0271             end
0272             %
0273             if isfield(launch,'epol_XYZ') && isfield(launch,'N_XYZ');
0274                 %
0275                 [Nperpp,Nperpm,Kxyz,ep_xyz,em_xyz] = colddisp_dke_jd(Npar0,launch.omega_rf,NaN,0,1,1,B0);
0276                 %
0277                 if isinf(equil.Rp),
0278                     BX0 = Bx0;
0279                     BY0 = -Bz0;
0280                 else
0281                     BX0 = Bx0*cos(z0/equil.Rp) - Bz0*sin(z0/equil.Rp);
0282                     BY0 = -Bx0*sin(z0/equil.Rp) - Bz0*cos(z0/equil.Rp);
0283                 end
0284                 BZ0 = By0;
0285                 %
0286                 B0_XYZ = [BX0;BY0;BZ0];
0287                 %
0288                 z_XYZ = B0_XYZ/B0;
0289                 y_XYZ = cross(B0_XYZ,launch.N_XYZ)/norm(cross(B0_XYZ,launch.N_XYZ));
0290                 x_XYZ = cross(y_XYZ,z_XYZ);
0291                 %
0292                 M_xyzXYZ = [x_XYZ,y_XYZ,z_XYZ];
0293                 ep_XYZ = M_xyzXYZ*ep_xyz;
0294                 em_XYZ = M_xyzXYZ*em_xyz;
0295                 %
0296                 rayinit.ap2 = abs(launch.epol_XYZ'*ep_XYZ).^2;
0297                 rayinit.am2 = abs(launch.epol_XYZ'*em_XYZ).^2;
0298             end
0299             %
0300         else
0301             rayinit.yrho0(iy) = NaN;%Initial radial position at launch
0302             rayinit.ytheta0(iy) = NaN;%Initial poloidal position at launch
0303             rayinit.yz0(iy) = NaN;%Initial axial position at launch
0304             rayinit.ym0(iy) = NaN;%Initial poloidal mode number
0305             rayinit.ykz0(iy) = NaN;%Initial toroidal mode number
0306             rayinit.yNpar0(iy) = NaN;%Initial index of refraction
0307             rayinit.ydNpar0(iy) = NaN;%initial Ray spectral width
0308             rayinit.yP0_2piRp(iy) = NaN;%Lineic initial power density initial power in the ray (W/m)
0309             rayinit.w0 = NaN;
0310             rayinit.z0 = NaN;
0311             rayinit.ap2 = NaN;
0312             rayinit.am2 = NaN;
0313             rayinit.yXfrac = NaN;
0314         end
0315     end
0316     %
0317 elseif strcmp(launch.type,'LH'),
0318     if isfield(launch,'bPlhtot') && ~isinf(equil.Rp),
0319         bP0_2piRp = launch.bPlhtot/(2*pi*equil.Rp);%initial power in the ray (W) [1,n_lobe]
0320     elseif isfield(launch,'bPlhtot_2piRp')
0321         bP0_2piRp = launch.bPlhtot_2piRp;
0322     else
0323         error('ERROR in main_rayinit_launch_jd.m -> inconsistent launch input: use bPlhtot if Rp is finite only, or bPlhtot_2piRp otherwise.');
0324     end
0325     %
0326     bNpar0 = launch.bNpar0;%Initial Npar value [1,n_lobe]
0327     bdNpar0 = launch.bdNpar0;%ray spectral width for Fokker-Planck calculations [1,n_lobe]
0328     ry0 = launch.rZ0 - equil.Zp;
0329     rby0 = ry0(:)*ones(1,length(bP0_2piRp));%[m_pol,n_lobe]
0330     %
0331     % Ray initial characteristics
0332     %
0333     if ~isfield(launch,'i_ref') || isempty(launch.i_ref),
0334         launch.i_ref = -1;%Initial radial position index, relative to LCFS [1,1]
0335     end
0336     %
0337     if launch.i_ref < 0,
0338         launch.i_ref = length(equil.psi_apRp) + launch.i_ref;%Initial radial position index [1,1]
0339     end
0340     %
0341     tmask = diff(equil.pty(launch.i_ref,:))*launch.LFS > 0;
0342     %
0343     rbx0 = interp1(equil.pty(launch.i_ref,tmask),equil.ptx(launch.i_ref,tmask),rby0,'spline');
0344     rbx0(rbx0 == 0) = eps;
0345     rbtheta0 = atan(rby0./rbx0) + pi*(rbx0 < 0) + 2*pi*(rbx0 > 0 & rby0 < 0);
0346     %
0347     rho0 = psi2rho_jd(equil,equil.psi_apRp(launch.i_ref)/equil.psi_apRp(end),'spline');%Initial radial position at launch%
0348     %
0349     rayinit.omega_rf = launch.omega_rf;%wave frequency
0350     %
0351     nb = size(rby0,2);% number of beams
0352     nr = size(rby0,1);% number of poloidal locations
0353     nm = length(launch.m0);% number of poloidal numbers
0354     %
0355     if all(imag(launch.m0) == 0),
0356         launch.m0 = launch.m0 + 1i/nm;%uniform power distribution
0357     end
0358     %
0359     % tail model
0360     %
0361     if ~isfield(launch,'tail'),% no tail model
0362         tail.mode = 0;
0363     else
0364         tail = launch.tail;
0365     end
0366     %
0367     for ib = 1:nb,% sum over lobes
0368         if imag(bNpar0(ib)) ~= 0,% case where Nz is prescribed instead of Npar - tail model only applies to Nz
0369             %
0370             bNy(ib) = real(bNpar0(ib));
0371             bNpar0(ib) = imag(bNpar0(ib));% nz
0372             %
0373         else
0374             bNy(ib) = NaN;
0375         end
0376     end
0377     %
0378     if tail.mode == 1,% Gaussian tail fluctuations
0379         %
0380         bfwhm = tail.bfwhm; 
0381         %
0382         if length(fwhm_npar0) == 1,% if tail bfwhm is the same for all lobes
0383             bfwhm = repmat(bfwhm,[1,nb]);
0384         end
0385         %
0386     elseif tail.mode == 2,% Linear tail fluctuations
0387         %
0388         bNparmax_tail = tail.bNparmax_tail;
0389         bopt_tail = real(tail.bopt_tail);
0390         %
0391         if length(bNparmax_tail) == 1,% if tail bNparmax is the same for all lobes
0392             bNparmax_tail = repmat(bNparmax_tail,[1,nb]);
0393         end
0394         %
0395         if length(bopt_tail) == 1,
0396             bopt_tail = repmat(bopt_tail,[1,nb]);
0397         end
0398         %
0399         bNparmax_tail(bNparmax_tail == 0) = 6.5/sqrt(equil.pTe(1));
0400         %
0401     elseif tail.mode == 3,% calc_tail_jd fluctuation or static model
0402         %
0403         bNparmax_tail = tail.bNparmax_tail;
0404         bn_tail = tail.bn_tail;
0405         bP_tail = tail.bP_tail;
0406         bopt_tail = tail.bopt_tail;
0407         %
0408         if length(bNparmax_tail) == 1,% if tail bNparmax is the same for all lobes
0409             bNparmax_tail = repmat(bNparmax_tail,[1,nb]);
0410         end
0411         %
0412         bNparmax_tail(bNparmax_tail == 0) = 6.5/sqrt(equil.pTe(1));
0413         %
0414         if length(bn_tail) == 1,
0415             bn_tail = repmat(bn_tail,[1,nb]);
0416         end
0417         %
0418         if length(bP_tail) == 1,
0419             bP_tail = repmat(bP_tail,[1,nb]);
0420         end
0421         %
0422         if length(bopt_tail) == 1,
0423             bopt_tail = repmat(bopt_tail,[1,nb]);
0424         end
0425         %
0426         for ib = 1:nb,% sum over lobes
0427             %
0428             [btNpar0_tail{ib},btdNpar0_tail{ib},btP0_2piRp_tail{ib}] = calc_tail_jd(bNpar0(ib),bdNpar0(ib),bP0_2piRp(ib),bNparmax_tail(ib),bn_tail(ib),bP_tail(ib),bopt_tail(ib));
0429             %
0430         end
0431         %
0432     end
0433     %
0434     iy = 0;
0435     %
0436     for ib = 1:nb,% sum over initial lobes
0437         for ir = 1:nr,% sum over vertical positions
0438             if tail.mode == 1,% Gaussian tail fluctuations
0439                 %
0440                 if ~isinf(tail.dtn),% fluctuating spectrum model
0441                     if bfwhm(ib) > 0,
0442                         tNpar0_tail = randg_yp(1,bNpar0(ib),bfwhm(ib));%New random npar0 around npar_antenna
0443                     else
0444                         tNpar0_tail = randg_yp(1,bNpar0(ib),bdNpar0(ib));%New random npar0 around npar_antenna
0445                     end   
0446                     %
0447                     tdNpar0_tail = bdNpar0(ib);
0448                     tP0_2piRp_tail = bP0_2piRp(ib);
0449                     %
0450                     disp(['Gaussian distribution -> npar_antenna: ',num2str(bNpar0(ib)),', new npar0: ',num2str(tNpar0_tail),' with Gaussian distribution FWHM = ',num2str(bfwhm(ib))]);
0451                 else
0452                     error('static model for gaussian probability distribution not yet implemented.')%% TODO
0453                 end
0454                 %
0455             elseif tail.mode == 2,% Linear tail fluctuations
0456                 %
0457                 if ~isinf(tail.dtn),% fluctuating spectrum model
0458                     if sign(bNpar0(ib)) > 0 && bNpar0(ib) < bNparmax_tail(ib),
0459                         tNpar0_tail = randlin_yp(1,bopt_tail(ib),[bNpar0(ib),bNparmax_tail(ib)]);%New random npar0
0460                     elseif sign(bNpar0(ib)) < 0 && bNpar0(ib) > - bNparmax_tail(ib),
0461                         tNpar0_tail = randlin_yp(1,bopt_tail(ib),[-bNparmax_tail(ib),bNpar0(ib)]);%New random npar0
0462                     end  
0463                     %
0464                     tdNpar0_tail = bdNpar0(ib);
0465                     tP0_2piRp_tail = bP0_2piRp(ib);
0466                     %
0467                     disp(['Linear distribution -> npar_antenna: ',num2str(bNpar0(ib)),', new npar0: ',num2str(tNpar0_tail),' with linear distribution h = ',num2str(bopt_tail(ib)),', npar_max = ',num2str(bNparmax_tail(ib))]);
0468                 else
0469                     error('Linear probability distribution only valid for fluctuating model.')
0470                 end
0471                 %
0472             elseif tail.mode == 3,% calc_tail_jd fluctuations or static model
0473                 %
0474                 if ~isinf(tail.dtn),% fluctuating spectrum model
0475                     %
0476                     if length(btNpar0_tail{ib}) > 1,
0477                         tNpar0_tail = randreject_yp(@fcalc_tail_yp,{btNpar0_tail{ib},btdNpar0_tail{ib},btP0_2piRp_tail{ib},bP0_2piRp(ib),[btNpar0_tail{ib}(1)-2*btdNpar0_tail{ib}(1)*sign(btNpar0_tail{ib}(1)),btNpar0_tail{ib}(end)+2.0*btdNpar0_tail{ib}(1)*sign(btNpar0_tail{ib}(end))]},1,0);
0478                     else
0479                         tNpar0_tail = randg_yp(1,bNpar0(ib),bdNpar0(ib));%New random npar0 around npar_antenna Gaussian shape)
0480                         %bNpar0_tail = bNpar0(ib);%npar0 = npar_antenna
0481                     end
0482                     %
0483                     tdNpar0_tail = bdNpar0(ib);
0484                     tP0_2piRp_tail = bP0_2piRp(ib);
0485                     %
0486                     disp(['Calculated distribution -> npar_antenna: ',num2str(bNpar0(ib)),', new npar0: ',num2str(tNpar0_tail),' with linear distribution h = ',num2str(bopt_tail(ib)),', npar_max = ',num2str(bNparmax_tail(ib))]);
0487                 else
0488                     % static spectrum model
0489                     tNpar0_tail = btNpar0_tail{ib};
0490                     tdNpar0_tail = btdNpar0_tail{ib};
0491                     tP0_2piRp_tail = btP0_2piRp_tail{ib};
0492                     %
0493                 end
0494             else% no tail model
0495                 tNpar0_tail = bNpar0(ib);
0496                 tdNpar0_tail = bdNpar0(ib);
0497                 tP0_2piRp_tail = bP0_2piRp(ib);                
0498             end
0499             %
0500             for itail = 1:length(tNpar0_tail),% sum of tail lobes
0501                 %
0502                 if ~isnan(bNy(ib)),% Npar0 is in fact Nz readjust Npar0 with Ny
0503                     tNpar0_tail(itail) = bNy(ib) + 1i*tNpar0_tail(itail);
0504                 end
0505                 %
0506                 for im = 1:nm,% sum over initial m
0507                     %
0508                     iy = iy + 1;
0509                     %
0510                     rayinit.yrho0(iy) = rho0;%Initial radial position at launch
0511                     rayinit.ytheta0(iy) = rbtheta0(ir,ib);%Initial poloidal position at launch
0512                     if isfield(launch,'z0')
0513                         rayinit.yz0(iy) = launch.z0;%Initial axial position at launch
0514                            rayinit.ykz0(iy) = launch.kz0;%Initial axial wave vector
0515                     elseif isfield(launch,'phi0') && ~isinf(equil.Rp),
0516                         rayinit.yphi0(iy) = launch.phi0;%Initial toroidal position at launch
0517                         rayinit.yn0(iy) = launch.n0;%Initial toroidal mode number
0518                     else
0519                         error('ERROR in main_rayinit_launch_jd.m -> inconsistent launch input: use phi0 if Rp is finite only, or z0 otherwise.');
0520                     end
0521                     %
0522                     rayinit.ym0(iy) = real(launch.m0(im));%Initial poloidal mode number
0523                     %
0524                     rayinit.yNpar0(iy) = tNpar0_tail(itail);%Initial index of refraction
0525                     rayinit.ydNpar0(iy) = tdNpar0_tail(itail);%initial Ray spectral width
0526                     rayinit.yP0_2piRp(iy) = tP0_2piRp_tail(itail)*imag(launch.m0(im))/nr;%Lineic initial power density initial power in the ray (W/m)
0527                     %
0528                     rayinit.yib(iy) = ib;
0529                     rayinit.yir(iy) = ir;
0530                     rayinit.yim(iy) = im;
0531                     rayinit.yit(iy) = itail;
0532                 end
0533             end
0534         end
0535     end
0536 else
0537     error('launch type not recognised')
0538 end
0539 %
0540 rayinit.launch = launch;
0541 %