optgrid_dke_jd

 This function generates a non-uniform radial grid of nr points based on
 the linear damping of waves such that there are more points in regions of
 large absorption.

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

0001 function [rho_S,pnmax_S] = optgrid_dke_jd(waves,rho_S,pnmax_S,equil,w_mask,display_mode)
0002 %
0003 % This function generates a non-uniform radial grid of nr points based on
0004 % the linear damping of waves such that there are more points in regions of
0005 % large absorption.
0006 %
0007 % by J. Decker (joan.decker@cea.fr) (CEA/DSM/IRFM) and Y. Peysson (yves.peysson@cea.fr) (CEA/DSM/IRFM)
0008 %
0009 if nargin < 6,
0010     display_mode = 0;
0011 end
0012 %
0013 if nargin < 5 || any(isnan(w_mask)),
0014     w_mask = 1:length(waves);
0015 elseif max(w_mask) > length(waves)
0016     info_dke_yp(2,'Warning : max(w_mask) > length(waves), dismissed.');
0017     iw = w_mask > length(waves);
0018     w_mask(iw) = [];
0019 end
0020 %
0021 if length(rho_S) ~= 1 || real(rho_S) <= 1,
0022     return
0023 else
0024     nr = real(rho_S);
0025     rho0 = imag(rho_S);
0026     clear rho_S
0027     info_dke_yp(2,'Radial grid automaticaly generated for optimizing the RF absorption.');
0028 end
0029 %
0030 nr_dep = 1000;%number of grid point in radial damping calculations
0031 finefac = 1;%fineness factor
0032 %
0033 dP = zeros(1,nr_dep);
0034 %
0035 rho_dep_S = linspace(0,1,nr_dep+1);
0036 drho_dep = 1/nr_dep;
0037 %
0038 for iw = w_mask;
0039     %
0040     if isfield(waves{iw},'rays'),
0041         for iray = 1:length(waves{iw}.rays),
0042             %
0043             ray = waves{iw}.rays{iray};
0044             %
0045             if isfield(ray,'sP_2piRp_lin')  && isfield(ray,'srho'),
0046                 if isfield(ray,'sP_2piRp_coll')
0047                     dP = dP + drho_dep*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin_coll,nr_dep);
0048                 else
0049                     dP = dP + drho_dep*radialdampingprofile_jd(ray.srho,ray.sP_2piRp_lin,nr_dep);
0050                 end
0051             end
0052             %
0053         end
0054     end
0055     %
0056 end
0057 %
0058 Ptot = sum(dP);
0059 %
0060 if Ptot == 0,
0061     dP0 = 1;
0062 else
0063     %
0064     if pnmax_S < 0,
0065         %
0066         pnmax_S = -pnmax_S;
0067         %
0068         rho_dep = (rho_dep_S(1:end-1) + rho_dep_S(2:end))/2;
0069         irmax = find(dP == max(dP));
0070         rhoabs = rho_dep(irmax(1));
0071         %
0072         prho = equil.ptx(:,1)/equil.ptx(end,1);
0073         pTe = equil.pTe;
0074         %
0075         te0 = max(pTe);
0076         teabs = interp1(prho,pTe,rhoabs);
0077         %
0078         pnmax_S = max([pnmax_S*sqrt(teabs/te0),5]);% 5 is minimum to avoid losses in maxwellian
0079     end
0080     %
0081     dP0 = Ptot/finefac/nr_dep;
0082 end
0083 %
0084 dP = dP + dP0*ones(1,nr_dep);
0085 rP = [0,cumsum(dP)];
0086 P = rP(end);
0087 %
0088 rho_S = interp1(rP,rho_dep_S,P*linspace(0,1,nr+1));
0089 %
0090 if rho0 > 0 && rho_S(2) > rho0,
0091     rho_S = [0,rho0,rho_S(2:end)];
0092 end
0093 %
0094 
0095 
0096 
0097 
0098 
0099 
0100 
0101 
0102 
0103 
0104 
0105 
0106 
0107