make_wave_TScirc_e1_LH_karney_offaxis_2

 script make_wave_TScirc_LHtest

 Parameters for test mode LHCD calculations
 This function has a benchmarking purpose only
 Approximations: 
   - Circular equilibrium required
   - The polarization is purely electrostatic: in that case, only the parallel polarization term is non-zero.
   - Small FLR limit: J0(z) -> 1. 

 by J. Decker (RLE/MIT) <jodecker@mit.edu> and Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr>

0001 % script make_wave_TScirc_LHtest
0002 %
0003 % Parameters for test mode LHCD calculations
0004 % This function has a benchmarking purpose only
0005 % Approximations:
0006 %   - Circular equilibrium required
0007 %   - The polarization is purely electrostatic: in that case, only the parallel polarization term is non-zero.
0008 %   - Small FLR limit: J0(z) -> 1.
0009 %
0010 % by J. Decker (RLE/MIT) <jodecker@mit.edu> and Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr>
0011 %
0012 clear all
0013 %
0014 id_wave = 'LH_karney_offaxis_2';
0015 %
0016 id_equil = 'TScirc_e1';%For plasma equilibrium
0017 path_equil = '../EQUIL_files/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0018 %
0019 [equil] = load_structures_yp('equil',id_equil,path_equil);
0020 %
0021 omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit
0022 %Option parameter for cross-comparison between old LH code:
0023 %    - (1): 1/vpar dependence
0024 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0025 opt_lh = 2; % [1,1]
0026 %
0027 %
0028 blist_lh = [1]; %Index of selected scenarios for the LH wave. [1,n_scenario_lh]
0029 %
0030 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0031 % or (Npar_lh,dNpar_lh) for Gaussian shape
0032 %
0033 norm_ref = 0;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0034 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0035 yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0036 yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0037 %
0038 yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0039 yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0040 yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0041 ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0042 %
0043 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0044 yD0_in_c_lh = [0.5];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0045 %
0046 yD0_in_lh_prof = [1];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0047 ypeak_lh = [0.40];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0048 ywidth_lh = [0.10];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0049 %
0050 ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0051 %
0052 % -------------------------------------------------------------------------
0053 %
0054 % parameters for makefile
0055 %
0056 ns = 1000;%for sufficient radial resolution (even number)
0057 method = 'spline';
0058 %
0059 %==========================================================================
0060 %
0061 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;
0062 %
0063 [equilDKE] = equilibrium_jd(equil);
0064 %
0065 srho = 1:-1/ns:1/ns;%radial position along the rays
0066 %
0067 Te0 = equilDKE.xTe(1);
0068 ne0 = equilDKE.xne(1);
0069 lnc_e_ref = 31.3 - 0.5*log(ne0) + log(Te0*1000);%Coulomb logarithm (Sauter et al. Phys. Plasmas, 6 (1999) 2834)
0070 betath_ref = sqrt(Te0/mc2);
0071 nhuth_ref = qe^4*ne0.*lnc_e_ref./(4*pi*e0^2*me^2.*(clum*betath_ref).^3);%Electron collision frequency (s-1)
0072 %
0073 sTe = interp1(equilDKE.xrho,equilDKE.xTe,srho,method);
0074 sne = interp1(equilDKE.xrho,equilDKE.xne,srho,method);
0075 slnc_e = 31.3 - 0.5*log(sne) + log(sTe*1000);%Coulomb logarithm (Sauter et al. Phys. Plasmas, 6 (1999) 2834)
0076 sbetath = sqrt(sTe/mc2);
0077 snhuth = qe^4*sne.*slnc_e./(4*pi*e0^2*me^2.*(clum*sbetath).^3);%Electron collision frequency (s-1)
0078 swpe = sqrt(qe^2*sne/me/e0);
0079 %
0080 if norm_ref == 0,
0081     srbetath = ones(1,ns);
0082     srnhuth = ones(1,ns);
0083 else
0084     srbetath = sbetath/betath_ref;
0085     srnhuth = snhuth/nhuth_ref;
0086 end
0087 %
0088 nb_lh = length(yD0_in_c_lh);% number of stored scenarios
0089 %
0090 sythetab_lh = ones(ns,1)*ythetab_lh;
0091 %
0092 syvparmin_lh = ones(ns,1)*yvparmin_lh;
0093 syvparmax_lh = ones(ns,1)*yvparmax_lh;
0094 syNparmin_lh = ones(ns,1)*yNparmin_lh;
0095 syNparmax_lh = ones(ns,1)*yNparmax_lh;
0096 syNpar_lh = ones(ns,1)*yNpar_lh;
0097 sydNpar_lh = ones(ns,1)*ydNpar_lh;
0098 syD0_in_c_lh = ones(ns,1)*yD0_in_c_lh;
0099 sypeak_lh = ones(ns,1)*ypeak_lh;
0100 sywidth_lh = ones(ns,1)*ywidth_lh;
0101 %
0102 syrbetath_lh = repmat(srbetath',[1,nb_lh]);
0103 syrnhuth_lh = repmat(srnhuth',[1,nb_lh]);
0104 sybetath_loc =  repmat(sbetath',[1,nb_lh]);
0105 %
0106 ytest_maskvpar_lh =  ~isnan(yvparmin_lh) & ~isnan(yvparmax_lh);
0107 ytest_maskNpars_lh =  ~isnan(yNparmin_lh) & ~isnan(yNparmax_lh);
0108 ytest_maskNparg_lh =  ~isnan(yNpar_lh) & ~isnan(ydNpar_lh);
0109 ytest_input_lh = [ytest_maskvpar_lh;ytest_maskNpars_lh;ytest_maskNparg_lh];
0110 %
0111 if sum(ytest_input_lh) ~= ones(1,nb_lh)
0112     error('LH Fourier spectrum input not compatible')
0113 end
0114 %
0115 syvparmin_lh_loc = syvparmin_lh./syrbetath_lh;%Local Normalization (when norm_ref = 1)
0116 syvparmax_lh_loc = syvparmax_lh./syrbetath_lh;%Local Normalization (when norm_ref = 1)
0117 %
0118 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)
0119 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)
0120 ytest_maskNpars_lh = logical(ytest_maskNpars_lh + ytest_maskvpar_lh);
0121 %
0122 syNpar_lh(:,ytest_maskNpars_lh) = (syNparmax_lh(:,ytest_maskvpar_lh) + syNparmin_lh(:,ytest_maskvpar_lh))/2;
0123 %
0124 % sydNpar_lh is made imaginary to identify square spectrum type (see rfwave)
0125 sydNpar_lh(:,ytest_maskNpars_lh) = i*(syNparmax_lh(:,ytest_maskvpar_lh) - syNparmin_lh(:,ytest_maskvpar_lh));
0126 %
0127 yD0_in_lh_prof = logical(yD0_in_lh_prof);
0128 %
0129 syrho = srho'*ones(1,nb_lh);
0130 %
0131 syD0_in_lh = zeros(ns,nb_lh);
0132 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
0133 syD0_in_lh(:,~yD0_in_lh_prof) = syD0_in_c_lh(:,~yD0_in_lh_prof);%Uniform radial profile
0134 %
0135 syD0_lh = syD0_in_lh./(syrbetath_lh.*syrbetath_lh.*syrnhuth_lh);%Local Normalization
0136 %
0137 % Approximate values for relevant polarization and normalized power flow
0138 % (Note: These eventually cancel out, but are useful for correct equivalent power scaling)
0139 %
0140 syepol_z_lh = 1./swpe'*omega_lh*ones(1,nb_lh);%approx Kperp=1, cancel out anyway
0141 syphi_z_lh = 1./syNpar_lh;%approx Kperp=1, cancel out anyway
0142 %
0143 syepol_p_lh = [0]*ones(ns,nb_lh);%Not relevant in small FLR limit
0144 syepol_m_lh = [0]*ones(ns,nb_lh);%Not relevant in small FLR limit
0145 syphi_x_lh = [0]*ones(ns,nb_lh);%Not relevant
0146 syphi_y_lh = [0]*ones(ns,nb_lh);%Not relevant
0147 %
0148 %Normalization of the diffusion coefficient for comparison with old calc. (Karney, Shoucri), see rfdiff_dke_jd
0149 %
0150 synorm_lh = (srho*equilDKE.ap*equilDKE.Rp*me.*slnc_e.*swpe.^2*omega_lh)'*ones(1,nb_lh).*abs(syphi_z_lh)./syepol_z_lh.^2;
0151 %
0152 syP0_lh = zeros(ns,nb_lh);
0153 syP0_lh(:,ytest_maskNpars_lh) = synorm_lh(:,ytest_maskNpars_lh).*syD0_lh(:,ytest_maskNpars_lh).*(syNparmax_lh(:,ytest_maskNpars_lh) - syNparmin_lh(:,ytest_maskNpars_lh));%Square power spectrum
0154 syP0_lh(:,ytest_maskNparg_lh) = synorm_lh(:,ytest_maskNparg_lh).*syD0_lh(:,ytest_maskNparg_lh).*sydNpar_lh(:,ytest_maskNparg_lh)*sqrt(pi);%Gaussien power spectrum
0155 %
0156 % Note: The initial power in the ray is a scalar. In order to take into
0157 % account the location dependence calculated above, and which results from
0158 % D0 normalization, profiles and dnpar, we include this dependence in phi,
0159 % while P0 is taken to be its maximum value along the ray.
0160 % (since Dql ~ P/|phi|)
0161 %
0162 yP0_lh = max(syP0_lh);
0163 syphi_z_lh = syphi_z_lh.*(ones(ns,1)*yP0_lh)./(syP0_lh+eps);
0164 %
0165 %
0166 syNperp_lh = syNpar_lh./syepol_z_lh;%Electrostatic limit. Indicative, not relevant in small FLR limit
0167 %
0168 syR = equilDKE.Rp + equilDKE.ap*syrho.*cos(sythetab_lh);
0169 syZ = equilDKE.Zp + equilDKE.ap*syrho.*sin(sythetab_lh);
0170 syphi = [0]*ones(ns,nb_lh);% so that the calculated finc is 1.
0171 %
0172 for i_lh = blist_lh,
0173     ray.P0_2piRp = yP0_lh(i_lh)'/(2*pi*equilDKE.Rp);% Power launched in the ray (W) [1,1]
0174         %
0175     ray.sx = syR(:,i_lh)' - equil.Rp;% Ray location (major radius) along the propagation [1,ns]
0176     ray.sphi = syZ(:,i_lh)';% Ray location (toroidal angle) along the propagation [1,ns]
0177     ray.sy = syphi(:,i_lh)' - equil.Zp;% Ray location (vertical) along the propagation [1,ns]
0178         %
0179     ray.ss = [];% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0180     ray.spsin = [];% Ray location (normalized magnetic flux) along the propagation [1,ns]
0181     ray.stheta = [];% Ray location (poloidal angle) along the propagation [1,ns]
0182         %
0183     ray.skx = [];% Ray wave vector (major radius) along the propagation [1,ns]
0184     ray.sn = [];% Ray wave number (toroidal) along the propagation [1,ns]
0185     ray.sky = [];% Ray wave vector (vertical) along the propagation [1,ns]
0186         %
0187     ray.sNpar = syNpar_lh(:,i_lh)';% Ray parallel wave number [1,ns]
0188     ray.sNperp = syNperp_lh(:,i_lh)';% Ray perpendicular wave number [1,ns]
0189     ray.sdNpar = reshape(sydNpar_lh(:,i_lh),[1,ns]);% Ray spectral width in Npar [1,ns]
0190         %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar)
0191         %
0192     syepol_pmz_lh = [syepol_p_lh(:,i_lh)';syepol_m_lh(:,i_lh)';syepol_z_lh(:,i_lh)'];
0193     syphi_xyz_lh = [syphi_x_lh(:,i_lh)';syphi_y_lh(:,i_lh)';syphi_z_lh(:,i_lh)'];
0194         %
0195     ray.sepol_pmz = syepol_pmz_lh;% polarization vector [3,ns] complex
0196     ray.sphi_xyz = syphi_xyz_lh;% normalized power flow [3,ns]
0197     %
0198     rays{i_lh} = ray;
0199 end
0200 %
0201 wave.omega_rf = omega_lh;
0202 wave.opt_rf = opt_lh;%option for FLR calculations [1,1]
0203 wave.opt_disp = 0;  % option for solving the wave equation, which gives (ray.sepol_pmz, ray.sphi, ray.sNperp)
0204 wave.equil_id = id_equil;
0205 wave.rays = rays;
0206 %
0207 save_str = ['WAVE_',id_equil,'_',id_wave,'.mat'];
0208 eval(['save ',save_str,' wave']);