make_beam_wave_jd

 This function generates a beam made of nd x nchi power tubes along a central ray

 Note : ray.swidth can be a scalar in which case the beam divergence is assumed to be zero.

0001 function wave = make_beam_wave_jd(wave,equil,C3POdisplay,nd,nchi,w0,z0,nsm)
0002 %
0003 % This function generates a beam made of nd x nchi power tubes along a central ray
0004 %
0005 % Note : ray.swidth can be a scalar in which case the beam divergence is assumed to be zero.
0006 %
0007 if nargin < 8,
0008     nsm = 1;%smoothing parameter
0009 end
0010 %
0011 if nargin < 7,
0012     z0 = NaN;
0013 end
0014 %
0015 if nargin < 6,
0016     w0 = NaN;
0017 end
0018 %
0019 if nargin < 5,
0020     nchi = 1;
0021 end
0022 %
0023 if nargin < 4,
0024     nd = 1;
0025 end
0026 %
0027 Nrays = length(wave.rays);
0028 %
0029 rays = cell(1,Nrays*nd*nchi);
0030 %
0031 for iray = 1:Nrays,
0032     %
0033     ray = wave.rays{iray};
0034     %
0035     ss = ray.ss;
0036     sphi_perp = abs(ray.sphi_xyz(1,:));
0037     sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0038     omega_rf = wave.omega_rf;
0039     %
0040     ns = length(ss);
0041     %
0042     clum = 299792458;%Speed of light (m/s)
0043     %
0044     if ~isnan(w0),% beam width information given
0045         if ~isnan(z0)% gaussian beam representation
0046             zr = omega_rf.*w0^2/clum/2;
0047             swidthfac = sqrt(1 + (ss - z0).^2/zr^2);
0048         else
0049             swidthfac = ones(1,ns);
0050         end
0051         %
0052         swidth = w0*swidthfac/sqrt(2);
0053         %
0054         %ray.sdNpar = clum*sphi_perp./swidth./sphi_norm/omega_rf;
0055         %
0056     else
0057         %warning('No beam width specified. Beam width determined from dNpar')
0058         swidth = clum*sphi_perp./ray.sdNpar./sphi_norm/omega_rf;
0059     end
0060     %
0061     if nd == 1 || nchi == 1,
0062         %
0063         rays{iray} = ray;
0064         %
0065         if ~isnan(swidth),
0066             rays{iray}.swidth = swidth;
0067         end
0068         %
0069         continue
0070     end
0071     %
0072     sx = ray.sx;
0073     sy = ray.sy;
0074     sphi = ray.sphi;
0075     sNpar = ray.sNpar; 
0076     sNperp = ray.sNperp;
0077     %
0078     if nsm > 1,
0079         sx = smooth(ss,sx,nsm).';
0080         sy = smooth(ss,sy,nsm).';
0081         sphi = smooth(ss,sphi,nsm).';
0082         %
0083         if length(sNpar) == ns,
0084             sNpar = smooth(ss,sNpar,nsm).';
0085         end
0086         if length(sNperp) == ns,
0087             sNperp = smooth(ss,sNperp,nsm).';
0088         end
0089     end
0090     %
0091     sR = sx + equil.Rp;
0092     %
0093     sdx = zeros(1,ns);
0094     sdy = zeros(1,ns);
0095     sdphi = zeros(1,ns);
0096     %
0097     sdx(2:ns-1) = sx(3:ns) - sx(1:ns-2);
0098     sdy(2:ns-1) = sy(3:ns) - sy(1:ns-2);
0099     sdphi(2:ns-1) = sphi(3:ns) - sphi(1:ns-2);
0100     %
0101     sdx(1) = 2*(sx(2) - sx(1));
0102     sdx(ns) = 2*(sx(ns) - sx(ns-1));
0103     %
0104     sdy(1) = 2*(sy(2) - sy(1));
0105     sdy(ns) = 2*(sy(ns) - sy(ns-1));
0106     %
0107     sdphi(1) = 2*(sphi(2) - sphi(1));
0108     sdphi(ns) = 2*(sphi(ns) - sphi(ns-1));
0109     %
0110     sv = sqrt(sdx.^2 + sR.^2.*sdphi.^2 + sdy.^2);
0111     sz = sqrt(sdx.^2 + sR.^2.*sdphi.^2);
0112     %
0113     for id = 1:nd,
0114         %
0115         sdi = swidth*sqrt(log(2*nd./(2*(nd - id) + 1)));
0116         %
0117         for jchi = 1:nchi,
0118             %
0119             ir = jchi + nchi*(id - 1 + (iray-1)*nd);
0120             %
0121             schij = (2*jchi - 1)*pi/nchi;
0122             %
0123             sxR_ij = -sdi.*(cos(schij).*sdx.*sdy./sv - sin(schij).*sR.*sdphi)./sz;
0124             sxy_ij = sdi.*cos(schij).*(sdx.^2 + sR.^2.*sdphi.^2)./sz./sv;
0125             sxphi_ij = -sdi.*(cos(schij).*sR.*sdphi.*sdy./sv + sin(schij).*sdx)./sz;
0126             %
0127             sR_ij = sqrt((sR + sxR_ij).^2 + sxphi_ij.^2);
0128             sy_ij = sy + sxy_ij;
0129             sphi_ij = sphi + atan(sxphi_ij./(sR + sxR_ij));
0130             %
0131             sx_ij = sR_ij - equil.Rp;
0132             %
0133             rays{ir}.P0_2piRp = ray.P0_2piRp/(nd*nchi);
0134             rays{ir}.sx = sx_ij;
0135             rays{ir}.sphi = sphi_ij;
0136             rays{ir}.sy = sy_ij;
0137             %
0138             rays{ir}.ss = [];
0139             rays{ir}.spsin = [];
0140             rays{ir}.stheta = [];
0141             %
0142             % The spectral properties must be the same as the central ray
0143             %
0144             rays{ir}.skx = ray.skx;
0145             rays{ir}.sn = ray.sn;
0146             rays{ir}.sky = ray.sky;
0147             rays{ir}.sm = ray.sm;
0148             rays{ir}.sn = ray.sn;
0149             rays{ir}.sNpar = sNpar;
0150             rays{ir}.sNperp= sNperp;
0151             rays{ir}.sdNpar = ray.sdNpar;
0152             rays{ir}.sepol_pmz = ray.sepol_pmz;
0153             rays{ir}.sphi_xyz = ray.sphi_xyz;
0154             %
0155             if ~isnan(swidth),
0156                 rays{ir}.swidth = swidth;
0157             end
0158             %
0159             rays{ir}.rayinits = ray.rayinits;
0160             %
0161         end
0162     end
0163 end
0164 %
0165 wave = rmfield(wave,'rays');
0166 wave.rays = rays;
0167 %
0168 if isfield(wave,'rayparam') && wave.rayparam.rel_opt == 0,
0169     wave.opt_disp = real(wave.opt_disp) - 1i;
0170 else
0171     wave.opt_disp = real(wave.opt_disp) + 1i;
0172 end    
0173 %
0174 
0175 
0176