NSTX support, JSON input file, Colored messages, bug fixes

D3D/d3d_equil.pro

@@ -9,13 +9,13 @@ FUNCTION d3d_equil,inputs,grid,det
     gfile=inputs.profile_dir+shot_str+'/g'+profile_str
 	gfiletest=findfile(gfile)
 	if gfiletest ne '' then begin
-		print,'RESTORING EQUILIBRIUM FROM GFILE'
+		print,'Restoring equilbrium from gfile'
 		g=readg(gfile) 
 	endif else begin
-		print,'FETCHING EQUILIBRIUM FROM MDS+'
+		print,'Fetching equilbrium from MDS+'
 		g=readg(inputs.shot,inputs.time*1000,RUNID=inputs.equil,status=gerr)
 		if gerr ne 1 then begin
-			print,'READG FAILED'
+			print,'readg failed'
 			equil={err:1}
 			goto,GET_OUT
 		endif

D3D/d3d_input.json

@@ -0,0 +1,104 @@
+//-----------------------------------------------------
+//                PREFIDA INPUT FILE
+//-----------------------------------------------------
+{
+"shot":146088,          // Shot Number
+"time":1.385,           // Time 
+"runid":"146088H05",    // runid of FIDASIM
+"device":"D3D",         // D3D,NSTX,AUGD,MAST
+"result_dir":"/u/stagnerl/FIDASIM/RESULTS/D3D/",           
+                        // Location where results will be stored
+"profile_dir":"/u/heidbrin/OANB/AUG/",                     
+                        // Location of profile files
+
+//-----------------------------------------------------
+// Fast-ion distribution function from transp
+//-----------------------------------------------------
+"cdf_file":"/e/alfven/FIDAsim/D3D/146088/146088H02_fi_9.cdf",    
+                        // CDF file from transp with the distribution funciton
+"emin":0.0,             // Minimum energy used from the distribution function
+"emax":100.0,           // Maximum energy used from the distribution function
+"pmin":-1.0,            // Minimum pitch used from the distribution function
+"pmax":1.0,             // Maximum pitch used from the distribution function
+
+//-----------------------------------------------------
+// Beam/FIDA/EQUILIBRIUM Selection
+//-----------------------------------------------------
+"isource":5,            // Beam source index (FIDASIM only simulates one NBI source)
+"einj":0.,              // [keV] If 0, get data from MDS+
+"pinj":0.,              // [MW] If 0, get data from MDS+
+"diag":"OBLIQUE",       // Name of the diagnostics 
+"equil":"EFIT01",       // Name of equilibrium. Ex. for D3D EFIT02
+
+//-----------------------------------------------------
+// Discharge Parameters
+//-----------------------------------------------------
+"btipsign":-1.0,        // Bt and Ip are in the opposite direction   
+"ab":2.01410178,        // Atomic mass of beam [u]
+"ai":2.01410178,        // Atomic mass of hydrogenic plasma ions [u]
+"impurity_charge":6,    // 5: BORON, 6: carbon, 7: Nitrogen
+
+//-----------------------------------------------------
+// Wavelength Grid
+//-----------------------------------------------------
+"lambdamin":647.0,      // Minimum wavelength of wavelength grid [nm] 
+"lambdamax":667.0,      // Maximum wavelength of wavelength grid [nm] 
+"nlambda":2000,         // Number of wavelengths
+"dlambda":0.01,         // Wavelength seperation
+
+//---------------------------------------------------
+// Define FIDASIM grid in machine coordinates(x,y,z)
+//---------------------------------------------------
+"nx":40,                // Number of cells in x direction
+"ny":60,                // Number of cells in y direction
+"nz":50,                // Number of cells in z direction
+"xmin":-170.0,          // Minimum x value
+"xmax":-70.0,           // Maximum x value
+"ymin":-195.0,          // Minimum y value
+"ymax":-80.0,           // Maximum y value
+"zmin":-70.0,           // Minimum z value
+"zmax":70.0,            // Maximum z value
+
+"origin":[0,0,0],       // If using different a coordinate system, this is the origin 
+                        // in machine coordinates of the new system
+
+"alpha":0.0,            // Rotation angle in radians from +x about z axis that transforms machine
+                        // coordinates to the new system. 
+"beta":0.0,             // Rotation about +y axis
+
+//--------------------------------------------------
+// Define number of Monte Carlo particles
+//--------------------------------------------------
+"nr_fast":5000000,      // FIDA
+"nr_nbi":50000,         // Beam emission
+"nr_halo":500000,       // Halo contribution
+
+//--------------------------------------------------
+// Calculation of the weight function
+//--------------------------------------------------
+"ne_wght":50,               // Number of Energies 
+"np_wght":50,               // Number of Pitches 
+"nphi_wght":50,             // Number of Gyro-angles 
+"emax_wght":125.0,          // Maximum energy (keV)
+"ichan_wght":-1,            // -1 for all channels, otherwise a given channel index
+"dwav_wght":0.2,            // Wavelength interval
+"wavel_start_wght":651.0,   // Minimum wavelength
+"wavel_end_wght":663.0,     // Maximum wavelength
+
+//-------------------------------------------------
+// Simulation switches
+//-------------------------------------------------
+"calc_npa":0,           // (0 or 1) If 1 do a simulation for NPA
+"calc_spec":1,          // (0 or 1) If 1 then spectra is calculated
+"calc_birth":1,         // (0 or 1) If 1 then the birth profile is calculated
+"calc_brems":0,         // (0 or 1) If 0 use the IDL bremstrahlung calculation
+"calc_fida_wght":1,     // (0 or 1) If 1 then fida weight functions are calculated
+"calc_npa_wght":0,      // (0 or 1) If 1 then npa weight functions are calculated
+"load_neutrals":0,      // (0 or 1) If 1 then the neutral density is loaded from an existing run 
+"load_fbm":1            // (0 or 1) If 1 then the fbm is loaded (calc_spec/npa overwrites)
+
+//------------------------------------------------
+// Extra Variables
+//------------------------------------------------
+
+}

D3D/d3d_input.pro

@@ -55,12 +55,12 @@ beta=0.0
 nx=40				;; Number of cells in x direction
 ny=60				;; Number of cells in y direction
 nz=50				;; Number of cells in z direction
-xdim1=-520.			;; Minimum x value
-xdim2=-420.			;; Maximum x value
-ydim1=-40.			;; Minimum y value
-ydim2=40.			;; Maximum y value
-zdim1=-40.			;; Minimum z value
-zdim2=40.			;; Maximum z value
+xmin=-520.			;; Minimum x value
+xmax=-420.			;; Maximum x value
+ymin=-40.			;; Minimum y value
+ymax=40.			;; Maximum y value
+zmin=-40.			;; Minimum z value
+zmax=40.			;; Maximum z value
 
 ;;--------------------------------------------------
 ;; Define number of Monte Carlo particles
@@ -101,7 +101,7 @@ inputs={shot:shot,time:time,runid:runid,device:strupcase(device),install_dir:ins
        isource:isource,einj:einj,pinj:pinj,diag:diag,equil:equil,$
        btipsign:btipsign,ab:ab,ai:ai,impurity_charge:impurity_charge,$
        lambdamin:lambdamin,lambdamax:lambdamax,nlambda:nlambda,dlambda:dlambda,origin:origin,alpha:alpha,beta:beta,$
-       nx:nx,ny:ny,nz:nz,xdim1:xdim1,xdim2:xdim2,ydim1:ydim1,ydim2:ydim2,zdim1:zdim1,zdim2:zdim2,$
+       nx:nx,ny:ny,nz:nz,xmin:xmin,xmax:xmax,ymin:ymin,ymax:ymax,zmin:zmin,zmax:zmax,$
        nr_fida:nr_fida,nr_ndmc:nr_ndmc,nr_halo:nr_halo,ne_wght:ne_wght,np_wght:np_wght,nphi_wght:nphi_wght,$
        emax_wght:emax_wght,ichan_wght:ichan_wght,dwav_wght:dwav_wght,wavel_start_wght:wavel_start_wght,$
 	   wavel_end_wght:wavel_end_wght,calc_npa:calc_npa,calc_spec:calc_spec,calc_birth:calc_birth, $

D3D/d3d_profiles.pro

@@ -17,7 +17,12 @@ FUNCTION d3d_profiles,inputs,save=save
 	shot_str=strtrim(string(inputs.shot),1)
 	profile_str=shot_str+'.'+time_str
 
-	dir=inputs.profile_dir+shot_str+'/'
+	dir = inputs.profiles_dir
+    ;; Look for electron density in this directory
+    test = FILE_SEARCH(STRJOIN([dir,'dne*']))
+    ;; if we don't find it, then look in the shot sub-directory
+    if not test THEN dir = STRJOIN([dir,shot_str])+'/'
+
 	ne_string=dir+'dne'+profile_str
 	te_string=dir+'dte'+profile_str
 	ti_string=dir+'dti'+profile_str

LIB/read_json.pro

@@ -0,0 +1,93 @@
+FUNCTION tokenize,str,regex
+	;;; Splits string into tokens according to regex ;;;
+    json=str
+    group=[]
+    start=[]
+    finish=[]
+    cnt=0
+    while stregex(json,regex,/BOOLEAN) do begin
+        pos=stregex(json,regex,length=len)
+        tmp=strmid(json,pos,len)
+        group=[group,tmp]
+        if cnt eq 0 then begin
+            start=[start,pos]
+            finish=[finish,pos+len]
+        endif else begin
+            start=[start,finish[cnt-1]+pos]
+            finish=[finish,finish[cnt-1]+pos+len]
+        endelse
+        json=strmid(json,pos+len,strlen(json))
+        cnt+=1
+    endwhile
+    return,{group:group,start:start,finish:finish}   
+END
+
+FUNCTION json_minify,str
+    ;;; Strips out C like comments out of JSON string ;;;
+    ;;; Based of JSON.minify https://github.com/getify/JSON.minify ;;;
+    regex='"|(/\*)|(\*/)|(//)|'+string(10b)+'|'+string(13b)
+    slashes='(\\)*$'
+
+    in_string = 0
+    in_multi = 0 
+    in_single = 0
+    new_str=''
+    index=0
+    tokenizer=tokenize(str,regex)
+
+    for i=0,n_elements(tokenizer.group)-1 do begin        
+        if not (in_multi or in_single) then begin
+             tmp=strmid(str,index,tokenizer.start[i]-index)
+             if not in_string then begin
+                 tmp=strcompress(tmp,/remove_all)
+                 tmp=strjoin(strsplit(tmp,string(10b),/extract,/regex),'')
+                 tmp=strjoin(strsplit(tmp,string(13b),/extract,/regex),'')                
+             endif
+             new_str+=tmp
+        endif
+
+        index=tokenizer.finish[i]
+        val=tokenizer.group[i]
+
+        if val eq '"' and not (in_multi or in_single) then begin
+            escaped = stregex(strmid(str,0,tokenizer.start[i]),slashes,/EXTRACT)
+            if not in_string or (escaped eq '' or strlen(escaped) mod 2 eq 0) then in_string = not in_string
+            index-=1
+
+        endif else if not (in_string or in_multi or in_single) then begin
+            if val eq '/*' then in_multi = 1 else if val eq '//' then in_single=1
+
+        endif else if val eq '*/' and in_multi and not (in_string or in_single) then begin
+            in_multi=0
+        endif else if (stregex(val,string(10b),/boolean) or stregex(val,string(13b),/boolean)) $
+                      and not (in_multi or in_string) and in_single then begin
+            in_single=0
+        endif else if not ((in_multi or in_single) or $
+                      (stregex(val,string(10b),/boolean) or $ 
+                       stregex(val,string(13b),/boolean) or $
+                       stregex(val,string(9b),/boolean)  or $
+                       stregex(val,string(32b),/boolean) )) then begin
+            new_str+=val
+        endif
+    endfor            
+    new_str+=strmid(str,index,strlen(str)-index)
+
+    return,new_str
+END
+
+FUNCTION read_json, file
+    ;;; Reads a JSON file into a string ;;;
+    ;;; Strips out C like comments ;;;
+    ;;; Returns structure containing JSON variables ;;;
+    openr,lun,file,/GET_LUN
+    json=''
+    tmp=''
+    while not EOF(lun) do begin
+        readf,lun,tmp
+        json+=tmp
+        json+=string(10b)
+    endwhile
+    free_lun,lun
+    stripped_json=json_minify(json)
+    return,json_parse(stripped_json,/toarray,/tostruct)
+END

LIB/transp_geometry.pro

@@ -0,0 +1,50 @@
+;;PROCEDURE: transp_geometry
+;;DESCRIPTION: calculates xyz_src and xyz_pos from TRANSP geometry variables
+;;INPUTS:
+;;  RTCENA: Radius of tangency point [cm]
+;;  XLBAPA: Distance from center of beam source grid to aperture [cm]
+;;  XLBTNA: Distance from center of beam source grid to tangency point
+;;  XBZETA: Torodial angle [deg] Positive angles defined to be in the counter-clockwise direction
+;;  XYBAPA: Elevation above/below vacuum vessel midplane of beam centerline at aperture [cm]
+;;  XYBSCA: Elevation above/below vacuum vessel midplane of center of beam source grid [cm]
+;;  NLCO:   1 for Co-beam, -1 for Counter-beam
+;;KEYWORDS:
+;;  angle: Angle to add to XBZETA to rotate the beams into correct coordinates Ex. 90 for D3D,145 for NSTX
+;;  plot: Plot the beam
+;;OUTPUTS:
+;;  xyz_src: Beam source position in machine coordinates [cm]
+;;  xyz_pos: Point on beam line in vessel, halfway between aperture and tangency point
+PRO transp_geometry,RTCENA,XLBAPA,XLBTNA,XBZETA,XYBAPA,XYBSCA,NLCO,xyz_src,xyz_pos,angle=angle,plot=plot
+
+    if not keyword_set(angle) then angle=0
+    if n_elements(RTCENA) gt 1 then begin
+        print,'TRANSP_GEOMETRY does not take array arguements'
+        goto,GET_OUT
+    endif
+
+    phi_s=(XBZETA+angle)*!DPI/180.0
+    zs=XYBSCA
+    za=XYBAPA
+    alpha=asin((zs-za)/XLBAPA)
+    pdst=XLBTNA*cos(alpha)
+    rs=sqrt(RTCENA^2.0 + pdst^2.0)
+    dat=XLBTNA-XLBAPA
+    pdat=dat*cos(alpha)
+    ra=sqrt(RTCENA^2.0 + pdat^2.0)
+    beta_s=acos(RTCENA/rs)
+    beta_a=acos(RTCENA/ra)
+    phi_a=phi_s+NLCO*(beta_s-beta_a)
+
+    transp_src=[rs*cos(phi_s),rs*sin(phi_s),zs]
+    xyz_src=[ra*cos(phi_a),ra*sin(phi_a),za]
+    dir=(xyz_src-transp_src)
+    dir=dir/sqrt(total(dir*dir))
+    xyz_pos=xyz_src+dir*(dat/2.0)
+
+    if keyword_set(plot) then begin
+        plot,[transp_src[0],xyz_pos[0]],[transp_src[1],xyz_pos[1]],$
+             xrange=1.05*[-xyz_src[0],xyz_src[0]],$
+             yrange=1.05*[-xyz_src[1],xyz_src[1]]
+    endif
+    GET_OUT:
+END