imake_wave_TORAY_jd

 This function makes the WAVE file from TORAY ray-tracing ascii outputs

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

0001 function wave = imake_wave_TORAY_jd(equil)
0002 %
0003 % This function makes the WAVE file from TORAY ray-tracing ascii outputs
0004 %
0005 % by J. Decker (RLE/MIT) <jodecker@mit.edu> and Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr>
0006 %
0007 wave = '';
0008 %
0009 % get data
0010 %
0011 if isfield(equil,'info') & isfield(equil.info,'shotnum') & isfield(equil.info,'shotime');
0012     ray_file = ['rays.toray_',equil.info.shotnum,'t',equil.info.shotime];
0013     geom_file = ['raygeom_',equil.info.shotnum,'t',equil.info.shotime];
0014 else
0015     ray_file = 'rays.toray_';
0016     geom_file = 'raygeom_';
0017 end                
0018 %
0019 [ray_file,ray_path] = uigetfile('rays.toray_*.*',['Please select the toray rays data file :'],ray_file);
0020 [geom_file,geom_path] = uigetfile('raygeom_*.*',['Please select the toray launching file :'],[ray_path,geom_file]);
0021 %
0022 if all(ray_file) == 0 | all(ray_file) == 0,
0023     return
0024 end
0025 %
0026 % [ny] = textread([geom_path,geom_file],'%n',1);%number of rays
0027 % [yX_L,yY_L,yZ_L,ybeta_L,yalpha_L,yfreq] = textread([geom_path,geom_file],'%n%n%n%n%n%n',ny+1);
0028 %
0029 fid = fopen([geom_path,geom_file]);
0030 ny = fscanf(fid,'%d',1);
0031 for iy = 1:ny,
0032     raydata = fscanf(fid,'%f',7);
0033     yX_L(iy) = raydata(1);
0034     yY_L(iy) = raydata(2);
0035     yZ_L(iy) = raydata(3);
0036     ybeta_L(iy) = raydata(4);
0037     yalpha_L(iy) = raydata(5);
0038     yfreq(iy) = raydata(6);
0039 end
0040 fclose(fid);
0041 %
0042 % group by beams
0043 %
0044 freq = -1;
0045 ib = 0;
0046 %
0047 for iy = 1:ny,
0048     if yfreq(iy) ~= freq | abs(beta_L - ybeta_L(iy)) > 10, %new beam
0049         ib = ib + 1;
0050         beam{ib} = iy;
0051         freq = yfreq(iy);
0052         beta_L = ybeta_L(iy);
0053     else, %same beam
0054         beam{ib} = [beam{ib},iy];
0055     end
0056     nb = ib;
0057 end
0058 %
0059 disp(['-----> The number of beams is ',num2str(nb)]);
0060 for ib = 1:nb,
0061     freq(ib) = yfreq(beam{ib}(1));
0062     beta_L(ib) = mean(ybeta_L(beam{ib}));
0063     disp(['-----> Beam # ',num2str(ib),', frequency : ',num2str(freq(ib)),', beta_L : ',num2str(beta_L(ib))]);
0064 end
0065 %
0066 % select the beam
0067 %
0068 ib = input('Select a beam number [1] : ');
0069 if isempty(ib),
0070     ib = 1;
0071 end
0072 if ~isnumeric(ib) | ~any(ib == 1:nb),
0073     disp('-----> Invalid beam selection');
0074     return
0075 end
0076 %
0077 freq_new = input(['Wave Frequency [',num2str(freq(ib)),']: ']);
0078 if isempty(freq_new) 
0079     freq_new = freq(ib);
0080 end
0081 if ~isnumeric(freq(ib)),
0082     disp('-----> Invalid frequency value');
0083     return
0084 end
0085 %
0086 % load ray-tracing data from TORAY
0087 %
0088 [sx,sy,sz,sPabs,ss,sNperp,sNpar] = textread([ray_path,ray_file],'%n%n%n%n%n%n%n');
0089 %
0090 % separate into individual rays
0091 %
0092 iy = 0;
0093 while ~isempty(ss),
0094     iy = iy + 1;
0095     isep = min(find(ss == -1));
0096     %
0097     syx{iy} = sx(1:isep-1);
0098     syy{iy} = sy(1:isep-1);
0099     syz{iy} = sz(1:isep-1);
0100     syPabs{iy} = sPabs(1:isep-1);
0101     sys{iy} = ss(1:isep-1);
0102     syNperp{iy} = sNperp(1:isep-1);
0103     syNpar{iy} = sNpar(1:isep-1);
0104     ns(iy) = isep-1;
0105     %
0106     % Remove double points
0107     %
0108     idb = find(diff(sys{iy}) == 0);
0109     %
0110     syx{iy}(idb) = [];
0111     syy{iy}(idb) = [];
0112     syz{iy}(idb) = [];
0113     syPabs{iy}(idb) = [];
0114     sys{iy}(idb) = [];
0115     syNperp{iy}(idb) = [];
0116     syNpar{iy}(idb) = [];
0117     ns(iy) = ns(iy) - length(idb);
0118     %
0119     sx(1:isep) = [];
0120     sy(1:isep) = [];
0121     sz(1:isep) = [];
0122     sPabs(1:isep) = [];
0123     ss(1:isep) = [];
0124     sNperp(1:isep) = [];
0125     sNpar(1:isep) = [];
0126     %
0127 end
0128 %
0129 if iy ~= ny,
0130     disp('-----> Ray propagation and geometry data have different numbers of rays.');
0131     return
0132 end
0133 %
0134 % filter rays for the beam
0135 %
0136 syx = syx(beam{ib});
0137 syy = syy(beam{ib});
0138 syz = syz(beam{ib});
0139 syPabs = syPabs(beam{ib});
0140 sys = sys(beam{ib});
0141 syNperp = syNperp(beam{ib});
0142 syNpar = syNpar(beam{ib});
0143 ns = ns(beam{ib});
0144 %
0145 ny = length(beam{ib});
0146 %
0147 % normalization
0148 %
0149 for iy = 1:ny,
0150     syx{iy} = syx{iy}/1e2;%cm -> m
0151     syy{iy} = syy{iy}/1e2;
0152     syz{iy} = syz{iy}/1e2;
0153     sys{iy} = sys{iy}/1e2;
0154     %
0155     syR{iy} = sqrt(syx{iy}.^2 + syy{iy}.^2);
0156     syZ{iy} = syz{iy};
0157     syphi{iy} = -atan(syy{iy}./syx{iy});%LUKE definition of phi
0158 end
0159 %
0160 % ray spectral width
0161 %
0162 dNpar0 = input('Ray spectral width [0.02] : ');
0163 if isempty(dNpar0),
0164     dNpar0 = 0.02;
0165 end
0166 if ~isnumeric(dNpar0),
0167     disp('-----> Invalid choice');
0168     return
0169 end
0170 %
0171 for iy = 1:ny,
0172     sydNpar{iy} = dNpar0*ones(ns(iy),1);
0173 end
0174 %
0175 % ray power
0176 %
0177 P0 = input('Total power in the beam [MW] : ');
0178 if isempty(P0) | ~isnumeric(P0),
0179     disp('-----> Invalid choice');
0180     return
0181 end
0182 %
0183 for iy = 1:ny,
0184     ray.P0_2piRp = 1e6*P0/ny/(2*pi*equil.Rp);% Power launched in the ray (W) [1,1]
0185     %
0186     ray.sx = syR{iy}.' - equil.Rp;% Ray location (major radius) along the propagation [1,ns]
0187     ray.sphi = syphi{iy}.';% Ray location (toroidal angle) along the propagation [1,ns]
0188     %ray.sy = syZ{iy}.' - equil.Zp;% Ray location (vertical) along the propagation [1,ns]
0189     ray.sy = syZ{iy}.';% toray data vs Z = Zp (!!!!! or Z = Zp_offset???)
0190     %
0191     ray.ss = sys{iy}.';% Ray propagation distance parameter for each ray (put to NaN once the ray is absorbed) [1,ns]
0192     ray.spsin = [];% Ray location (normalized magnetic flux) along the propagation [1,ns]
0193     ray.stheta = [];% Ray location (poloidal angle) along the propagation [1,ns]
0194     %
0195     ray.skx = [];% Ray wave vector (major radius) along the propagation [1,ns]
0196     ray.sn = [];% Ray wave number (toroidal) along the propagation [1,ns]
0197     ray.sky = [];% Ray wave vector (vertical) along the propagation [1,ns]
0198     %
0199     ray.sNpar = syNpar{iy}.';% Ray parallel wave number [1,ns]
0200     ray.sNperp = syNperp{iy}.';% Ray perpendicular wave number [1,ns]
0201     ray.sdNpar = sydNpar{iy}.';% Ray spectral width in Npar [1,ns]
0202     %       if sdNpar is real -> gaussian spectrum, if sdNpar is imag. -> square spectrum (centered on sNpar)
0203     %
0204     ray.sepol_pmz = [];% polarization vector [3,ns] complex
0205     ray.sphi_xyz = [];% normalized power flow [3,ns]
0206     %
0207     ray.sP_2piRp_lin = ray.P0_2piRp*(1 - syPabs{iy}.');
0208     %
0209     rays{iy} = ray;
0210 end
0211 %
0212 kmode = input('Plasma model for wave propagation : (0:cold, 1:warm, 2:kinetic nonrelativistic, 3:kinetic relativistic) [0] :');
0213 if isempty(kmode),
0214     kmode = 0;
0215 end
0216 if ~isnumeric(kmode) | ~any(kmode == [0:4]),
0217     disp('-----> Invalid plasma model option');
0218     return
0219 end
0220 %
0221 opt_rf = input('FLR effects : (0:all, 1:small) [0] :');
0222 if isempty(opt_rf),
0223     opt_rf = 0;
0224 end
0225 if ~isnumeric(opt_rf) | ~any(opt_rf == [0:1]),
0226     disp('-----> Invalid FLR effects option');
0227     return
0228 end
0229 %
0230 wave_id_o = ['TORAY_',num2str(ib)];
0231 %
0232 wave_id = input(['Please provide an id for the wave [',wave_id_o,'] : '],'s');
0233 if isempty(wave_id),
0234     wave_id = wave_id_o;
0235 end
0236 %
0237 wave.id = wave_id;
0238 wave.omega_rf = 2*pi*freq_new*1e9;
0239 wave.opt_rf = opt_rf;%option for FLR calculations [1,1]
0240 wave.opt_disp = 1 + kmode;  % option for solving the wave equation, which gives (ray.sepol_pmz, ray.sphi, ray.sNperp)
0241 wave.equil_id = equil.id;
0242 wave.kextra = 0;
0243 wave.rays = rays;
0244 %
0245 wave.info.type = 'RT';
0246 wave.info.RT = 'TORAY';
0247 wave.info.geom = [geom_path,geom_file];
0248 wave.info.rays = [ray_path,ray_file];
0249 %
0250 %