make_rayinit_genray_jd

 script make_rayinit_genray

 Creates rayinit structure from GENRAY output files

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

0001 % script make_rayinit_genray
0002 %
0003 % Creates rayinit structure from GENRAY output files
0004 %
0005 % by J. Decker (RLE/MIT) <jodecker@mit.edu> and Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr>
0006 %
0007 clear all
0008 close all
0009 clc
0010 %
0011 %
0012 id_rayinit = 'GENRAY';%
0013 %
0014 % Source file parameters
0015 %
0016 filename = 'genray.nc';
0017 path_file = '../GENRAY/Source_files/';
0018 %
0019 % Equilibrium parameters
0020 %
0021 id_equil = 'CMOD_060728014_01100';%For plasma equilibrium
0022 path_equil = '../../EQUIL/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0023 %
0024 % Wave parameters
0025 %
0026 m0 = NaN;
0027 %
0028 method = 'spline';
0029 %
0030 %==========================================================================
0031 %
0032 [qe,me,mp,mn,e0,mu0,re,mc2,clum] = pc_dke_yp;
0033 %
0034 % Load equilibrium
0035 %
0036 if exist([path_equil,'EQUIL_',id_equil,'.mat']) == 0,%Magnetic equilibrium
0037     info_dke_yp(2,[path_equil,'EQUIL_',id_equil,'.mat does not exist']);
0038     break
0039 else
0040     eval(['load ',path_equil,'EQUIL_',id_equil,'.mat']);
0041     info_dke_yp(2,'Structure ''equil'' successfully loaded');
0042 end
0043 %
0044 % load structure with NETCDF data
0045 %
0046 data_out_netcdf = load_netcdf_jd([path_file,filename]);
0047 %
0048 [id_equil_out,omega_rf,sys,sypsin,syR,syZ,syphi,...
0049           syNpar,sydNpar,syNperp,...
0050           syepol_x,syepol_y,syepol_z,syphi_x,...
0051           syalpha,syB,syne,syTe,...
0052           yP0_in,yR_in,yZ_in,yphi_in,yNtor_in,yNpol_in] = process_genray_jd(data_out_netcdf);
0053 %
0054 ny = length(yR_in);   
0055 %
0056 Rp = equil.Rp;
0057 Zp = equil.Zp;
0058 %
0059 yx_in = yR_in - Rp;
0060 yy_in = yZ_in - Zp;
0061 %
0062 yr_in = sqrt(yx_in.^2 + yy_in.^2);
0063 yx_in(yx_in == 0) = eps;%-eps;
0064 ytheta_in = atan(yy_in./yx_in) + pi*(yx_in < 0) + 2*pi*(yx_in > 0 & yy_in < 0);     
0065 %
0066 ptr = sqrt(equil.ptx.^2 + equil.pty.^2);
0067 %
0068 psi_apRp = equil.psi_apRp;
0069 psin = psi_apRp/psi_apRp(end);
0070 theta = equil.theta;
0071 %
0072 % interpolation for r at ray poloidal locations
0073 %
0074 pyr_in = interp1(theta,ptr.',ytheta_in,method).';
0075 %
0076 % interpolation for psi at ray position (along cst theta lines)
0077 %
0078 for iy = 1:ny,
0079     ypsin_in(iy) = interp1(pyr_in(:,iy),psin,yr_in(iy),method).';
0080 end
0081 %
0082 yrho_in = psi2rho_jd(equil,ypsin_in);
0083 %
0084 if isnan(m0),%calculate m0 from yNpol and n0 from yNtor
0085     yBx_in = interp2(psin,theta,equil.ptBx.',ypsin_in,ytheta_in,method);
0086     yBy_in = interp2(psin,theta,equil.ptBy.',ypsin_in,ytheta_in,method);
0087     yBPHI_in = interp2(psin,theta,equil.ptBPHI.',ypsin_in,ytheta_in,method);
0088     %
0089     yBP_in = sqrt(yBx_in.^2 + yBy_in.^2);
0090     yBT_in = yBPHI_in;
0091     yB_in = sqrt(yBP_in.^2 + yBT_in.^2);
0092     %
0093     yP_in = yBP_in./yB_in;
0094     yT_in = yBT_in./yB_in;
0095     %
0096     yNphi_in = yNtor_in;
0097     yn0_in = yR_in.*yNphi_in*omega_rf/clum;
0098     %
0099     ycalpha_in = (yBy_in.*yx_in - yBx_in.*yy_in)./yr_in./yBP_in;
0100     %
0101     if 0,%determine yNpar_in from ym0_in and yn0_in
0102         yNs_in = yNpol_in;
0103         ym0_in = yr_in.*yNs_in*omega_rf/clum./ycalpha_in;
0104         %
0105         yNpar_in = (ycalpha_in.*yP_in.*ym0_in./yr_in + sign(yBPHI_in).*yT_in.*yn0_in./yR_in)*clum/omega_rf;
0106         %yNpar_in = yNpol_in.*yP_in + yNtor_in.*sign(yBPHI_in).*yT_in;
0107         figure(1),clf
0108         plot(1:ny,yNpar_in,'b',1:ny,syNpar(1,:),'r--');
0109         legend({'From m0 and no','From syNpar(1,:)'});
0110     else %determine ym0_i from yNpar_inn and yn0_in
0111         yNpar_in = syNpar(1,:);
0112         %
0113         ym0_in = (yNpar_in*omega_rf/clum - sign(yBPHI_in).*yT_in.*yn0_in./yR_in).*yr_in./ycalpha_in./yP_in;
0114         %ym0_in = yr_in.*yNpol_in*omega_rf/clum./ycalpha_in;
0115         figure(1),clf
0116         plot(1:ny,yr_in.*yNpol_in*omega_rf/clum./ycalpha_in,'b',1:ny,ym0_in,'r--');
0117         legend({'From yNpol','From syNpar(1,:) and yNtor'});
0118     end        
0119     %
0120 else
0121     ym0_in = m0*ones(1,nb);
0122     yNpar_in = syNpar(1,:);
0123 end
0124 %
0125 % Identify the different beams
0126 %
0127 iy_peak = 1 + find(yP0_in(3:ny) <= yP0_in(2:ny-1) & yP0_in(1:ny-2) < yP0_in(2:ny-1));
0128 iy_dip = 1 + find(yP0_in(3:ny) >= yP0_in(2:ny-1) & yP0_in(1:ny-2) > yP0_in(2:ny-1));
0129 %
0130 nb = length(iy_peak);
0131 if length(iy_dip) ~= nb-1,
0132     error('individual peaks not determined')
0133 end
0134 %
0135 brho_in = yrho_in(iy_peak);
0136 btheta_in = ytheta_in(iy_peak);
0137 bphi_in = yphi_in(iy_peak);
0138 ym0_in = ym0_in(iy_peak);
0139 bNpar_in = yNpar_in(iy_peak);
0140 %
0141 bymin = [1,iy_dip];
0142 bymax = [iy_dip+1,ny];
0143 %
0144 for ib = 1:nb,
0145     bP0_in(ib) = sum(yP0_in(bymin(ib):bymax(ib)));
0146     ymin_dNpar = bymin(ib) - 1 + min(find(yP0_in(bymin(ib):bymax(ib)) >= yP0_in(iy_peak(ib))/exp(1)));
0147     ymax_dNpar = bymin(ib) - 1 + max(find(yP0_in(bymin(ib):bymax(ib)) >= yP0_in(iy_peak(ib))/exp(1)));
0148     bdNpar_in(ib) = (yNpar_in(ymax_dNpar + 1) - yNpar_in(ymin_dNpar - 1))/2;
0149 end
0150 %
0151 bP0_2piRp_in = bP0_in/(2*pi*Rp);
0152 %
0153 rayinit.omega_rf = omega_rf;
0154 rayinit.yrho0 = brho_in;
0155 rayinit.ytheta0 = btheta_in;
0156 rayinit.yphi0 = bphi_in;
0157 rayinit.ym0 = ym0_in;
0158 rayinit.yNpar0 = bNpar_in;
0159 rayinit.ydNpar0 = bdNpar_in;
0160 rayinit.yP0_2piRp = bP0_2piRp_in;
0161 %
0162 save_str = ['RAYINIT_',id_equil,'_',id_rayinit,'.mat'];
0163 eval(['save ',save_str,' rayinit']);
0164 
0165