;<pre><i>  example file for numerical data cubes

	; accessed in FORWARD by setting 'numcube'

	;If this is your first time using for_drive, then read
        ;<a href = "http://www.hao.ucar.edu/FORWARD/FOR_SSW/docs/INSTRUCTIONS/examples0_setup.html">Setup instructions</a> 

	; Your numerical cube should be already in the form expected by the FORWARD
	; codes, see 
        ;<a href = "http://www.hao.ucar.edu/FORWARD/FOR_SSW/idl/MODELS/NUMCUBE/make_my_cube.pro">FORWARD/MODELS_NUMCUBE/make_my_cube.pro</a> </i>

	; NOTE FOR CoMP CALCULATIONS, YOU MUST USE FULL PATH IN FILENAME
	; BECAUSE DIRECTORIES ARE CHANGED DURING COMP ANALYSIS
	; this is set below by keyword DIRNAME

	; if you don't have a cube yet, you can run this example with
	; default simulation "yfan2dcube" (set below).  Ordinarily, if
	; for_drive,'numcube' is run with no CUBENAME keyword set, 
	; there is a default cube1.dat thatlives in $FORWARD/MODELS/NUMCUBE 
	; directory. It is an axisymmetric model with a powerlaw falloff 
	; in density where B=Br=1 everywhere. 

	; default is to have no density or magnetic field outside of your
	; cube boundaries. You can change this with the keywords listed in
        ;<a href = "http://www.hao.ucar.edu/FORWARD/FOR_SSW/idl/MODELS/NUMCUBE/numcubeprams.pro">FORWARD/MODELS_NUMCUBE/numcubeprams.pro</a> </i>

	; Note on speed of numcube calculations: If you are running multiple
	; plots on the same cube a given SSWIDL session, you can run the
	; first for_drive as usual or with the keyword nreinit=1. You can then
	; use nreinit=0 in subsequent calls to the same cube. This will
	; make use of internal common blocks so that you don't have to load
	; your cube each time, it cuts down on computing time. NOTE: that
	; because of the common blocks, the numcube model can be quite fussy
	; so if you want to plot a different cube, restart SSWIDL or it
	; will probably puke.

	; Finally, this program makes some example plots: the integrated 
	; observables (pB, CoMP) only need to be run once and saved to a save 
	; set.  So, if you want to run this several times (or variants of it),
	; set keyword DOINT=0 and it will look for savesets for the 
	; integrated observables

	;
	; whether or not to do polarization plane-of-sky slices</i>
	;

default,runslice,0

	;
	;<i> whether or not to run comp calculations</i>
	;

default,runcomp,0

	;
	; <i>whether or not to do the comp integration
	; only do once to create save files
	; after that turn off for speed</i>
	;

default,doint,1

	;
	; <i>where the data cube is</i>
	;

default,dirname,'$FORWARD_DB/'
default,simname,'yfan2dcube' 
cuben=dirname+simname+'.dat'

	;
	; <i>these are the horizontal and vertical axis limits on the plots</i>
	;

default,xxmin,0.8
default,xxmax,1.7
default,yymin,-0.7
default,yymax,0.7

	;
	; <i>this is the resolution of the calculation</i>
	;

default,ngr,50

	;
	; <i>this is the screen size of the plots</i>
	;	

nwinx=1024
nwiny=1024

	; <i>this is the number of arrows (ngr/narrow is number of points to skip)</i>

narrow=30

if doint ne 0 then begin

	;<i> generate LOS integral for pB</i>

 for_drive, 'numcube', cubename=cuben, ngrid=ngr,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='pb',/wl,/savemap,mapname=simname+'_wl_fullint',/noplots

	; <i>generate Plane of sky slice for CoMP</i>

 if runcomp eq 1 then begin
  for_drive, 'numcube', cubename=cuben, ngrid=ngr,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='LoI',/comp,/pos,/verbose,/savemap,mapname=simname+'_comp_pos',/noplots

	; <i>generate LOS for CoMP</i>

  for_drive, 'numcube', cubename=cuben, ngrid=ngr,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='LoI',/comp,/verbose,/savemap,mapname=simname+'_comp_fullint',/noplots

 endif
endif

	;
	; <i>make plots</i>
	;

loadct,0

!p.multi=[0,2,2]

for_drive, 'numcube', cubename=cuben, ngrid=ngr,nwinx=nwinx,nwiny=nwiny,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='dens',/plotlog,usecolor=0,/field,bscale=-.2,narr=narrow,winnum=3

for_drive, 'numcube', cubename=cuben, /wl,line='pB',/plotlog,usecolor=0,/field,bscale=-.2,narr=narrow,/verbose,readmap=simname+'_wl_fullint',/noerase,nreinit=0

for_drive, 'numcube', cubename=cuben, ngrid=ngr,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='bx',usecolor=41,/field,bscale=-.2,narr=narrow,/noerase,nreinit=0

for_drive, 'numcube', cubename=cuben, ngrid=ngr,xxmin=xxmin,xxmax=xxmax,yymin=yymin,yymax=yymax,line='bpos',usecolor=0,/field,bscale=-.2,narr=narrow,/noerase,nreinit=0

if runcomp eq 1 then begin
for_drive, 'numcube', cubename=cuben, /comp,line='LoI',usecolor=0,readmap=simname+'_comp_pos',docont=1,winnum=2,nwinx=nwinx,nwiny=nwiny,/plotlog,title='LoI_los'

for_drive, 'numcube', cubename=cuben, /comp,line='VoI',nreinit=0,/fieldlines,/stklines,bscale=-.2,pscale=-.2,narr=narrow,readmap=simname+'_comp_pos',/noerase,docont=1,usecolor=41,quantmap=quantmap_v_poi,title='VoI_los'

for_drive, 'numcube', cubename=cuben, /comp,line='LoI',usecolor=0,readmap=simname+'_comp_fullint',/noerase,docont=1,/plotlog,title='LoI integrated'

for_drive, 'numcube', cubename=cuben, /comp,line='VoI',nreinit=0,/fieldlines,/stklines,bscale=-.2,pscale=-.2,narr=narrow,readmap=simname+'_comp_fullint',/noerase,docont=1,usecolor=41,quantmap=quantmap_v_int,title='VoI integrated'

endif

!p.multi=0

	; <i>Calculates the Stokes P/I for a series of thin slices along the line
	; of sight and save the maps. Bear in mind the default is to surround
	; simulation with zero plasma. To put in a spherically symmetric
	; isothermal background, we add the keywords hydro=1,te=1d6</i>

if doint ne 0 and runcomp ne 0 then begin
 for ixoffset=1,13 do $
 for_drive,'numcube',cubename=cuben,$
 ngrid=ngr,nwinx=512,nwiny=512,$
 xoffset=-.875+float(ixoffset)*0.125,$
 extratitle='X-offset='+strtrim(float(-.875+float(ixoffset)*0.125),2),$
 /savemap,mapname=simname+'_I_X-offset='+strtrim(-1.75+float(ixoffset)*0.25,2),$
 line='LoI',/tif,/fieldlines,bscale=-1,narrow=10,/verbose,/comp,/pos,$
 /plotlog,xxmin=.8,xxmax=2,yymin=-.7,yymax=.7,$
 hydro=1,te=1d6,bout=2,imin=-2,imax=2,/noplot
endif


	;<i>this takes the save files you made above and plots P/I with
	;fieldlines and stokeslines as well. It saves the plots as tifs.</i>

if runslice ne 0 and runcomp ne 0 then begin
 for ixoffset=1,13 do $
 for_drive,'numcube',cubename=cuben,$
 ngrid=ngr,nwinx=512,nwiny=512,$
 xoffset=-.875+float(ixoffset)*0.125,$
 extratitle='X-offset='+strtrim(float(-.875+float(ixoffset)*0.125),2),$
 readmap=simname+'_I_X-offset='+strtrim(-1.75+float(ixoffset)*0.25,2),$
 line='LoI',/tif,/fieldlines,bscale=-1,narrow=10,/verbose,/comp,/pos,$
 /plotlog,xxmin=.8,xxmax=2,yymin=-.7,yymax=.7,$
 hydro=1,te=1d6,bout=2,imin=-2,imax=2
endif
;</pre>
end