This file contains the procedure to read SPIRAL fast-ion distribution file
Reads SPIRAL guiding center fast-ion distribution file
file: SPIRAL output file
time: Time [s]
ntotal: Total number of fast-ions
e_range: Energy range of particles to consider
p_range: Pitch range of particles to consider
particle_weight: Set particle/marker weight such that sum(particle_weights) = ntotal: Defaults to ntotal
/nparticles
btipsign: Sign of the dot product between the current and magnetic field (Required)
Distribution structure
IDL> dist = read_spiral("./mc_159243H06_9",ntotal=1e19,btipsign=-1)
FUNCTION read_spiral_header, file openr,lun,file,/get_lun npart = 0 weight = 0.0 gc = 0 err = 0 nheader = 0 line = '' header = '' readf, lun, line while strmid(line,0,1) eq ';' do begin nheader = nheader + 1 header = header + line + string(10B) if stregex(line,"particle weight",/fold_case) ne -1 then begin weight = float(stregex(line,"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?",/extract)) endif if stregex(line,"number",/fold_case) ne -1 then begin npart = long(stregex(line,"[0-9]+",/extract)) endif if stregex(line,"guiding center",/fold_case) ne -1 then begin gc = 1 endif if stregex(line,"R_[m]",/fold_case) ne -1 then begin if n_elements(strsplit(line)) ne 6 then begin err = 1 endif endif line = '' readf, lun, line endwhile free_lun,lun if gc eq 0 then begin error,'Not a guiding center distribution',/halt endif if npart eq 0 then begin error,'Number of particles le 0',/halt endif return, {nheader:nheader, npart:npart, header:header,weight:weight} END FUNCTION finite_struct, s n = N_TAGS(s) for i=0,n-1 do begin if finite(s.(i)) ne 1 then return, 0 endfor return, 1 END ;+ This file contains the procedure to read SPIRAL fast-ion distribution file FUNCTION read_spiral,file, time=time, ntotal=ntotal, e_range=e_range, p_range=p_range, $ particle_weight=particle_weight, btipsign=btipsign ;+#read_spiral ;+Reads SPIRAL guiding center fast-ion distribution file ;+*** ;+##Arguments ;+ **file**: SPIRAL output file ;+ ;+##Keyword Arguments ;+ **time**: Time [s] ;+ ;+ **ntotal**: Total number of fast-ions ;+ ;+ **e_range**: Energy range of particles to consider ;+ ;+ **p_range**: Pitch range of particles to consider ;+ ;+ **particle_weight**: Set particle/marker weight such that sum(particle_weights) = ntotal: Defaults to `ntotal`/nparticles ;+ ;+ **btipsign**: Sign of the dot product between the current and magnetic field (Required) ;+ ;+##Return Value ;+Distribution structure ;+ ;+##Example Usage ;+```idl ;+IDL> dist = read_spiral("./mc_159243H06_9",ntotal=1e19,btipsign=-1) ;+``` if not keyword_set(btipsign) then begin error,'btipsign is not set.',/halt endif if not file_test(file) then begin error, 'Nonexistent file: '+file,/halt endif if not keyword_set(time) then begin print, 'WARNING: Time is not set for SPIRAL distribution. Setting to 0.0 [s]' time = 0.d0 endif header = read_spiral_header(file) npart = header.npart nhead = header.nheader pweight = header.weight openr,unit,file, /get_lun ; Read header head = strarr(nhead) readf, unit, head ; Read in data r = dblarr(npart) z = dblarr(npart) phi = dblarr(npart) energy = dblarr(npart) pitch = dblarr(npart) cnt = 0L for i=0,npart-1 do begin s = {r:double(0),phi:double(0),z:double(0),energy:double(0),pitch:double(0)} readf,unit,s if not finite_struct(s) then begin continue endif else begin r[cnt] = s.r*100 phi[cnt] = s.phi z[cnt] = s.z*100 energy[cnt] = s.energy pitch[cnt] = s.pitch*btipsign ; SPIRAL pitch is defined relative to current cnt = cnt+1 endelse endfor free_lun, unit npart = cnt r = r[0:npart-1] z = z[0:npart-1] phi = phi[0:npart-1] energy = energy[0:npart-1] pitch = pitch[0:npart-1] orbit_class = replicate(1,npart) if pweight ne 0.0 then begin weight = replicate(pweight,npart) endif else begin if not keyword_set(particle_weight) then begin if not keyword_set(ntotal) then begin print, 'WARNING: ntotal is not set. Setting arbitrarily to 1e19' ntotal = 1.d19 endif weight = replicate(Ntotal/float(npart), npart) endif else begin weight = replicate(particle_weight, npart) endelse endelse if not keyword_set(e_range) then begin e_range = [min(energy),max(energy)] endif if not keyword_set(p_range) then begin p_range = [min(pitch),max(pitch)] endif ww = where(energy ge e_range[0] and energy le e_range[1],nw) if nw eq 0 then begin error,'No particles fall in requested energy range',/halt endif wwp = where(pitch[ww] ge p_range[0] and pitch[ww] le p_range[1],nwp) if nwp eq 0 then begin error,'No particles fall in the requested pitch range',/halt endif ww = ww[wwp] nw = n_elements(ww) print,'Number of markers: ',npart print,'Number of markers in phase space: ',nw dist_struct = {type:2,time:time,data_source:file_expand_path(file), $ nparticle:long(nw),nclass:1,r:r[ww],z:z[ww],phi:phi[ww],$ energy:energy[ww],pitch:pitch[ww],class:orbit_class[ww],$ weight:weight[ww]} return, dist_struct END