wave_process_jd

 This function completes the missing wave characteristics. In particular, it calculates the wave absorption coefficient and the wave linear damping.

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

0001 function wave = wave_process_jd(equil,wave,n_rf_list,method)
0002 %
0003 % This function completes the missing wave characteristics. In particular, it calculates the wave absorption coefficient and the wave linear damping.
0004 %
0005 % by J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr> and Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr>
0006 %
0007 if ~isfield(wave,'rays'),% Not a RT code (FW, etc)
0008     return
0009 end
0010 %
0011 if nargin < 4,
0012     method = 'linear';%'spline';%interpolation technique in psi
0013 end
0014 %
0015 % Equilibrium data with Fourier-enhanced poloidal precision and original radial grid
0016 %
0017 [equilDKE_Fourier] = equilibrium_jd(equil);
0018 %
0019 % Absorption calculation including dnpar effects
0020 %
0021 if ~any(isnan(n_rf_list)),
0022     if isfield(wave,'rayparam') && wave.rayparam.rel_opt == 0,
0023         wave.opt_disp = real(wave.opt_disp) - 1i;
0024     else
0025         wave.opt_disp = real(wave.opt_disp) + 1i;
0026     end
0027 end
0028 %
0029 if isfield(wave,'a_sdNpar') && wave.a_sdNpar > 0,
0030     %
0031     n_rf_list = 0;% used for Landau damping only
0032     %
0033     if isfield(wave,'rayparam') && wave.rayparam.rel_opt == 0,
0034         wave.opt_disp = real(wave.opt_disp) - (1 + wave.a_sdNpar)*1i;
0035     else
0036         wave.opt_disp = real(wave.opt_disp) + (1 + wave.a_sdNpar)*1i;
0037     end
0038 end
0039 %
0040 % Ray processing
0041 %
0042 for ib = 1:length(wave.rays),
0043     %
0044     if isfield(wave,'colldamp'),
0045         colldamp = wave.colldamp ;
0046     else
0047         if isfield(wave.rayparam,'colldamp'),%backward compatibility
0048             colldamp = wave.rayparam.colldamp;
0049         else
0050             colldamp = 0;%no collisional damping
0051         end
0052         %
0053         wave.colldamp = colldamp;
0054     end
0055     %
0056     ray = rayprocess_jd(wave.rays{ib},equilDKE_Fourier,wave.opt_disp,wave.omega_rf,n_rf_list,method,colldamp);
0057     %
0058     % Ray linear damping
0059     %
0060     salphaphi_lin = ray.salphaphi_lin;
0061     if any(~isnan(salphaphi_lin)),%absorption coefficient was calculated
0062         %
0063         sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0064         salpha_lin = salphaphi_lin./sphi_norm;
0065         salpha_lin(isnan(salpha_lin)) = 0;
0066         %
0067         if ~isfield(ray,'stau_lin') || all(isnan(ray.stau_lin)) || abs(imag(wave.opt_disp)) >= 1,
0068             %
0069             salpha_lin_h = (salpha_lin(1:end-1) + salpha_lin(2:end))/2;
0070             sds = diff(ray.ss);
0071             ray.stau_lin = [0,cumsum(salpha_lin_h.*sds)];
0072             %
0073         end
0074         %
0075         ray.sP_2piRp_lin = ray.P0_2piRp*exp(-ray.stau_lin);
0076         ray.sdP_2piRp_ds_lin = salpha_lin.*ray.sP_2piRp_lin;
0077         %
0078     elseif isfield(ray,'sP_2piRp_lin'),%power along the ray was calculated
0079         %
0080         stau_lin = -log(ray.sP_2piRp_lin/ray.sP_2piRp_lin(1));
0081         %
0082         ss_h = (ray.ss(1:end - 1) + ray.ss(2:end))/2;
0083         sdtauds_h_lin = diff(stau_lin)./diff(ray.ss);
0084         salpha_lin = interp1(ss_h,sdtauds_h_lin,ray.ss,'linear','extrap');
0085         %
0086         sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0087         salphaphi_lin = salpha_lin.*sphi_norm;
0088         salphaphi_lin(isnan(salphaphi_lin)) = 0;
0089         %
0090         ray.stau_lin = stau_lin;
0091         ray.salphaphi_lin = salphaphi_lin;
0092         ray.sdP_2piRp_ds_lin = salpha_lin.*ray.sP_2piRp_lin;
0093     end
0094     %
0095     % Collisional ray damping
0096     %
0097     if isfield(ray,'salphaphi_coll'),
0098         %
0099         salphaphi_coll = ray.salphaphi_coll;
0100         %
0101         if any(~isnan(salphaphi_coll)),%collisional damping
0102             %
0103             sphi_norm = sqrt(abs(ray.sphi_xyz(1,:)).^2 + abs(ray.sphi_xyz(2,:)).^2 + abs(ray.sphi_xyz(3,:)).^2);
0104             salpha_coll = salphaphi_coll./sphi_norm;
0105             salpha_coll(isnan(salpha_coll)) = 0;
0106             %
0107             if ~isfield(ray,'stau_coll') || all(isnan(ray.stau_coll)) || abs(imag(wave.opt_disp)) >= 1,
0108                 %
0109                 salpha_coll_h = (salpha_coll(1:end-1) + salpha_coll(2:end))/2;
0110                 sds = diff(ray.ss);
0111                 ray.stau_coll = [0,cumsum(salpha_coll_h.*sds)];
0112                 %
0113             end
0114             %
0115             ray.sP_2piRp_coll = ray.P0_2piRp*exp(-ray.stau_coll);%Power lost for collisional damping alone
0116             ray.sdP_2piRp_ds_coll = salpha_coll.*ray.sP_2piRp_coll;%Power lost for collisional damping alone
0117             %
0118         end
0119     end
0120     %
0121     if isfield(ray,'stau_lin'),
0122         %
0123         if isfield(ray,'stau_coll')
0124             ray.sP_2piRp_lin_coll = ray.P0_2piRp*exp(-ray.stau_lin-ray.stau_coll);%Combined linear + collisional damping
0125             ray.sdP_2piRp_ds_lin_coll = (salpha_lin + salpha_coll).*ray.sP_2piRp_lin_coll;%Combined linear + collisional damping
0126         else
0127             ray.sP_2piRp_lin_coll = ray.P0_2piRp*exp(-ray.stau_lin);%linear damping
0128             ray.sdP_2piRp_ds_lin_coll = salpha_lin.*ray.sP_2piRp_lin;%linear damping
0129         end
0130         %
0131     end
0132     %
0133     if isfield(ray,'stau_lin') && isfield(wave,'cut_fac'),
0134         %
0135         ns = length(ray.ss);
0136         %
0137         if isfield(ray,'stau_coll'),
0138             is_max = min([ns,wave.kextra+find(ray.stau_lin + ray.stau_coll > 10*wave.cut_fac)]);
0139             if ray.stau_lin(end) + ray.stau_coll(end)  < 10*wave.cut_fac && ~isinf(wave.cut_fac),
0140                 disp('WARNING: The ray might be too short. Increase tfinal')
0141             end
0142         else
0143             is_max = min([ns,wave.kextra+find(ray.stau_lin > 10*wave.cut_fac)]);
0144             if ray.stau_lin(end) < 10*wave.cut_fac && ~isinf(wave.cut_fac),
0145                 disp('WARNING: The ray might be too short. Increase tfinal')
0146             end            
0147         end
0148         %
0149         frays = fieldnames(ray);
0150         for ifield = 1:length(frays),
0151             ns_field = length(ray.(frays{ifield}));
0152             if ns_field == ns,
0153                 ray.(frays{ifield}) = ray.(frays{ifield})(:,1:is_max);
0154             end
0155         end
0156     end
0157     %
0158     wave.rays{ib} = ray;
0159     info_dke_yp(2,['Damping of ray #',num2str(ib),' calculated']);    
0160 end
0161 %
0162 wave.opt_disp = 0;
0163 %