make_idealLHwave_jd

 Creates a wave structure based on ideal LHCD parameters

0001 function wave = make_idealLHwave_jd(equil,wavestruct)
0002 %
0003 % Creates a wave structure based on ideal LHCD parameters
0004 %
0005 %
0006 omega_lh = wavestruct.omega_lh;
0007 opt_lh = wavestruct.opt_lh;
0008 norm_ref = wavestruct.norm_ref;
0009 yvparmin_lh = wavestruct.yvparmin_lh;
0010 yvparmax_lh = wavestruct.yvparmax_lh;
0011 yNparmin_lh = wavestruct.yNparmin_lh;
0012 yNparmax_lh = wavestruct.yNparmax_lh;
0013 yNpar_lh = wavestruct.yNpar_lh;
0014 ydNpar_lh = wavestruct.ydNpar_lh;
0015 yD0_in_c_lh = wavestruct.yD0_in_c_lh;
0016 yD0_in_lh_prof = wavestruct.yD0_in_lh_prof;
0017 ypeak_lh = wavestruct.ypeak_lh;
0018 ywidth_lh = wavestruct.ywidth_lh;
0019 ythetab_lh = wavestruct.ythetab_lh;
0020 %
0021 % parameters for makefile
0022 %
0023 ns = 1001;%for sufficient radial resolution (even number)
0024 method = 'spline';
0025 %
0026 %==========================================================================
0027 %
0028 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;
0029 %
0030 [equilDKE] = equilibrium_jd(equil);
0031 %
0032 srho = linspace(1,0,ns);%radial position along the rays
0033 %
0034 Te0 = equilDKE.xTe(1);
0035 ne0 = equilDKE.xne(1);
0036 lnc_e_ref = 31.3 - 0.5*log(ne0) + log(Te0*1000);%Coulomb logarithm (Sauter et al. Phys. Plasmas, 6 (1999) 2834)
0037 betath_ref = sqrt(Te0/mc2);
0038 nhuth_ref = qe^4*ne0.*lnc_e_ref./(4*pi*e0^2*me^2.*(clum*betath_ref).^3);%Electron collision frequency (s-1)
0039 %
0040 sTe = interp1(equilDKE.xrho,equilDKE.xTe,srho,method);
0041 sne = interp1(equilDKE.xrho,equilDKE.xne,srho,method);
0042 slnc_e = 31.3 - 0.5*log(sne) + log(sTe*1000);%Coulomb logarithm (Sauter et al. Phys. Plasmas, 6 (1999) 2834)
0043 sbetath = sqrt(sTe/mc2);
0044 snhuth = qe^4*sne.*slnc_e./(4*pi*e0^2*me^2.*(clum*sbetath).^3);%Electron collision frequency (s-1)
0045 swpe = sqrt(qe^2*sne/me/e0);
0046 %
0047 if norm_ref == 0,
0048     srbetath = ones(1,ns);
0049     srnhuth = ones(1,ns);
0050 else
0051     srbetath = sbetath/betath_ref;
0052     srnhuth = snhuth/nhuth_ref;
0053 end
0054 %
0055 nb_lh = length(yD0_in_c_lh);% number of stored scenarios
0056 %
0057 sythetab_lh = ones(ns,1)*ythetab_lh;
0058 %
0059 syvparmin_lh = ones(ns,1)*yvparmin_lh;
0060 syvparmax_lh = ones(ns,1)*yvparmax_lh;
0061 syNparmin_lh = ones(ns,1)*yNparmin_lh;
0062 syNparmax_lh = ones(ns,1)*yNparmax_lh;
0063 syNpar_lh = ones(ns,1)*yNpar_lh;
0064 sydNpar_lh = ones(ns,1)*ydNpar_lh;
0065 syD0_in_c_lh = ones(ns,1)*yD0_in_c_lh;
0066 sypeak_lh = ones(ns,1)*ypeak_lh;
0067 sywidth_lh = ones(ns,1)*ywidth_lh;
0068 %
0069 syrbetath_lh = repmat(srbetath',[1,nb_lh]);
0070 syrnhuth_lh = repmat(srnhuth',[1,nb_lh]);
0071 sybetath_loc =  repmat(sbetath',[1,nb_lh]);
0072 %
0073 ytest_maskvpar_lh =  ~isnan(yvparmin_lh) & ~isnan(yvparmax_lh);
0074 ytest_maskNpars_lh =  ~isnan(yNparmin_lh) & ~isnan(yNparmax_lh);
0075 ytest_maskNparg_lh =  ~isnan(yNpar_lh) & ~isnan(ydNpar_lh);
0076 ytest_input_lh = [ytest_maskvpar_lh;ytest_maskNpars_lh;ytest_maskNparg_lh];
0077 %
0078 if sum(ytest_input_lh) ~= ones(1,nb_lh)
0079     error('LH Fourier spectrum input not compatible')
0080 end
0081 %
0082 syvparmin_lh_loc = syvparmin_lh./syrbetath_lh;%Local Normalization (when norm_ref = 1)
0083 syvparmax_lh_loc = syvparmax_lh./syrbetath_lh;%Local Normalization (when norm_ref = 1)
0084 %
0085 syNparmin_lh(:,ytest_maskvpar_lh) = 1./syvparmax_lh_loc(:,ytest_maskvpar_lh)./sybetath_loc(:,ytest_maskvpar_lh);%Conversion from vpar representation to Npar one (if yvparmax_lh not a NaN)
0086 syNparmax_lh(:,ytest_maskvpar_lh) = 1./syvparmin_lh_loc(:,ytest_maskvpar_lh)./sybetath_loc(:,ytest_maskvpar_lh);%Conversion from vpar representation to Npar one (if yvparmin_lh not a NaN)
0087 ytest_maskNpars_lh = logical(ytest_maskNpars_lh + ytest_maskvpar_lh);
0088 %
0089 syNpar_lh(:,ytest_maskNpars_lh) = (syNparmax_lh(:,ytest_maskvpar_lh) + syNparmin_lh(:,ytest_maskvpar_lh))/2;
0090 %
0091 % sydNpar_lh is made imaginary to identify square spectrum type (see rfwave)
0092 sydNpar_lh(:,ytest_maskNpars_lh) = i*(syNparmax_lh(:,ytest_maskvpar_lh) - syNparmin_lh(:,ytest_maskvpar_lh));
0093 %
0094 yD0_in_lh_prof = logical(yD0_in_lh_prof);
0095 %
0096 syrho = srho'*ones(1,nb_lh);
0097 %
0098 syD0_in_lh = zeros(ns,nb_lh);
0099 syD0_in_lh(:,yD0_in_lh_prof) = 2*sqrt(log(2))*exp(-(syrho(:,yD0_in_lh_prof) - sypeak_lh(:,yD0_in_lh_prof)).^2./(sywidth_lh(:,yD0_in_lh_prof)/2/sqrt(log(2))).^2).*syD0_in_c_lh(:,yD0_in_lh_prof)/sqrt(pi);%Gaussian radial profile on midplane
0100 syD0_in_lh(:,~yD0_in_lh_prof) = syD0_in_c_lh(:,~yD0_in_lh_prof);%Uniform radial profile
0101 %
0102 syD0_lh = syD0_in_lh./(syrbetath_lh.*syrbetath_lh.*syrnhuth_lh);%Local Normalization
0103 %
0104 % Approximate values for relevant polarization and normalized power flow
0105 % (Note: These eventually cancel out, but are useful for correct equivalent power scaling)
0106 %
0107 syepol_z_lh = 1./swpe'*omega_lh*ones(1,nb_lh);%approx Kperp=1, cancel out anyway
0108 syphi_z_lh = 1./syNpar_lh;%approx Kperp=1, cancel out anyway
0109 %
0110 syepol_p_lh = [0]*ones(ns,nb_lh);%Not relevant in small FLR limit
0111 syepol_m_lh = [0]*ones(ns,nb_lh);%Not relevant in small FLR limit
0112 syphi_x_lh = [0]*ones(ns,nb_lh);%Not relevant
0113 syphi_y_lh = [0]*ones(ns,nb_lh);%Not relevant
0114 %
0115 %Normalization of the diffusion coefficient for comparison with old calc. (Karney, Shoucri), see rfdiff_dke_jd
0116 %
0117 synorm_2piRp_lh = (srho*equilDKE.ap*me.*slnc_e.*swpe.^2*omega_lh)'*ones(1,nb_lh).*abs(syphi_z_lh)./syepol_z_lh.^2/2/pi;
0118 %
0119 syP0_2piRp_lh = zeros(ns,nb_lh);
0120 syP0_2piRp_lh(:,ytest_maskNpars_lh) = synorm_2piRp_lh(:,ytest_maskNpars_lh).*syD0_lh(:,ytest_maskNpars_lh).*(syNparmax_lh(:,ytest_maskNpars_lh) - syNparmin_lh(:,ytest_maskNpars_lh));%Square power spectrum
0121 syP0_2piRp_lh(:,ytest_maskNparg_lh) = synorm_2piRp_lh(:,ytest_maskNparg_lh).*syD0_lh(:,ytest_maskNparg_lh).*sydNpar_lh(:,ytest_maskNparg_lh)*sqrt(pi);%Gaussien power spectrum
0122 %
0123 % Note: The initial power in the ray is a scalar. In order to take into
0124 % account the location dependence calculated above, and which results from
0125 % D0 normalization, profiles and dnpar, we include this dependence in phi,
0126 % while P0 is taken to be its maximum value along the ray.
0127 % (since Dql ~ P/|phi|)
0128 %
0129 yP0_2piRp_lh = max(syP0_2piRp_lh);
0130 syphi_z_lh = syphi_z_lh.*(ones(ns,1)*yP0_2piRp_lh)./(syP0_2piRp_lh+eps);
0131 %
0132 %
0133 syNperp_lh = syNpar_lh./syepol_z_lh;%Electrostatic limit. Indicative, not relevant in small FLR limit
0134 %
0135 syx = equilDKE.ap*syrho.*cos(sythetab_lh);
0136 syy = equilDKE.ap*syrho.*sin(sythetab_lh);
0137 syz = [0]*ones(ns,nb_lh);% so that the calculated finc is 1.
0138 syphi = [0]*ones(ns,nb_lh);% so that the calculated finc is 1.
0139 %
0140 for i_lh = 1:length(ythetab_lh),
0141     ray.P0_2piRp = yP0_2piRp_lh(i_lh)';% Power launched in the ray (W) [1,1]
0142         %
0143     ray.sx = syx(:,i_lh)';% Ray location (major radius) along the propagation [1,ns]
0144     ray.sy = syy(:,i_lh)';% Ray location (toroidal angle) along the propagation [1,ns]
0145     ray.sz = syz(:,i_lh)';% Ray location (vertical) along the propagation [1,ns]
0146     ray.sphi = syphi(:,i_lh)';% Ray location (vertical) along the propagation [1,ns]
0147         %
0148     ray.ss = [];% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0149     ray.spsin = [];% Ray location (normalized magnetic flux) along the propagation [1,ns]
0150     ray.stheta = [];% Ray location (poloidal angle) along the propagation [1,ns]
0151         %
0152     ray.skx = [];% Ray wave vector (major radius) along the propagation [1,ns]
0153     ray.sky = [];% Ray wave vector (vertical) along the propagation [1,ns]
0154     ray.skz = [];% Ray wave number (toroidal) along the propagation [1,ns]
0155     ray.sn = [];% Ray wave number (toroidal) along the propagation [1,ns]
0156         %
0157     ray.sNpar = syNpar_lh(:,i_lh)';% Ray parallel wave number [1,ns]
0158     ray.sNperp = syNperp_lh(:,i_lh)';% Ray perpendicular wave number [1,ns]
0159     ray.sdNpar = reshape(sydNpar_lh(:,i_lh),[1,ns]);% Ray spectral width in Npar [1,ns]
0160         %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar)
0161         %
0162     syepol_pmz_lh = [syepol_p_lh(:,i_lh)';syepol_m_lh(:,i_lh)';syepol_z_lh(:,i_lh)'];
0163     syphi_xyz_lh = [syphi_x_lh(:,i_lh)';syphi_y_lh(:,i_lh)';syphi_z_lh(:,i_lh)'];
0164         %
0165     ray.sepol_pmz = syepol_pmz_lh;% polarization vector [3,ns] complex
0166     ray.sphi_xyz = syphi_xyz_lh;% normalized power flow [3,ns]
0167     %
0168     rays{i_lh} = ray;
0169 end
0170 %
0171 wave.id = 'ideal';
0172 wave.omega_rf = omega_lh;
0173 wave.opt_rf = opt_lh;%option for FLR calculations [1,1]
0174 wave.opt_disp = 0;  % option for solving the wave equation, which gives (ray.sepol_pmz, ray.sphi, ray.sNperp)
0175 wave.equil_id = equil.id;
0176 wave.rays = rays;