Reads guiding center Monte Carlo NUBEAM fast-ion distribution file
infile: NUBEAM Monte Carlo distribution file
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("./spiral_file.TXT",time=1.0, ntotal=1e19)
FUNCTION read_mc_nubeam,infile,ntotal=ntotal,e_range=e_range,p_range=p_range,particle_weight = particle_weight,btipsign=btipsign ;+#read_mc_nubeam ;+Reads guiding center Monte Carlo NUBEAM fast-ion distribution file ;+*** ;+##Arguments ;+ **infile**: NUBEAM Monte Carlo distribution file ;+ ;+##Keyword Arguments ;+ **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("./spiral_file.TXT",time=1.0, ntotal=1e19) ;+``` if not keyword_set(btipsign) then begin error,'btipsign is not set.',/halt endif if not keyword_set(ntotal) and not keyword_set(particle_weight) then begin warn,'Ntotal not set. Setting arbitrarily to 1e19' ntotal = 1.d19 endif zzz=FINDFILE(infile) if zzz eq '' then begin error, 'Nonexistent file: '+infile,/halt endif openr,unit,infile, /get_lun line=' ' readf,unit,line ; read string readf,unit,line pos=strpos(line,'N=') if pos eq -1 then begin error,'Second line is missing the number of points',/halt endif w = stregex(line,'N *= *([0-9]*)',/sub,/extract) npts = long(w[1]) if npts lt 5 then begin error,'Too few points '+strtrim(npts,2),/halt endif ; Get time readf,unit,line parts=str_sep(line, ' ') w=where(parts eq 'TIME',nw) if nw eq 0 then begin error,'Time not found on 3rd line',/halt endif i=1 while 1 do begin s = parts[w[0]+i] if s ne '' and s ne '=' then begin time = float(s) break endif i = i+1 endwhile data=fltarr(4,npts) ready=0 while not ready do begin readf,unit,line ; Description line pos=strpos(line,'R(cm)') if pos gt -1 then ready=1 endwhile for i=0L,npts-1 do begin readf,unit,line ; read string line=strcompress(line) parts = str_sep(line, ' ') & np=n_elements(parts) if parts(0) eq '' then parts=shift(parts,-1) ; accommodate blank before 1st while parts(np-1) eq '' do np=np-1 if np ne 4 then begin error,'Wrong number of entries on line: '+line endif else begin data[*,i]=parts[0:np-1] endelse print,format='(f7.2,"%",A,$)',100.0*(i+1)/float(npts),string(13b) endfor free_lun, unit r=double(reform(data[0,0:npts-1])) w=double(reform(data[1,0:npts-1])) pitch=double(reform(data[2,0:npts-1]))*btipsign energy=reform(data[3,0:npts-1])*1.d-3 ;keV orbit_class = replicate(1,npts) if not keyword_set(particle_weight) then begin particle_weight = ntotal/float(npts) endif weight = replicate(particle_weight,npts) 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, 'Time: ',time print,'Number of markers: ',npts print,'Number of markers in phase space: ',nw print, 'Total Number of Fast-ions in phase space: ',particle_weight*nw fbm_struct = {type:2,time:double(time),data_source:file_expand_path(infile), $ nparticle:long(nw),nclass:1,r:r[ww],z:w[ww],$ energy:energy[ww],pitch:pitch[ww],class:orbit_class[ww],$ weight:weight[ww]} return, fbm_struct END