PRO flight_sim, my_obs_id, RUN_MARX=run_marx, $
	NUMRAYS=numrays, $
	POST_PROCESS=post_process, N_POST_EVENTS=n_post_events, IMAGE=image, $
	NEW_SCATTER = new_scatter, $
	SKIP_SETUP=skip_setup, MARX_SIM_DIR=marx_sim_dir
;
; Procedure to do a MARX simulation of a Flight Observation
;
; INPUTS:
;     my_obs_id:  'O-HAS-Capella-1.001'
;
;                O = on-orbit
;
; KEYWORDS:
;   RUN_MARX:
;   NUMRAYS:  value for the MARX NumRays parameter, default 100000.
;   POST_PROCESS:
;   N_POST_EVENTS: Returns the number of post-prcessed events in
;                  the idlsav file (so calling routine can re-run
;                  to get desired number of counts...)
;   IMAGE:  Make a cute color image...
;   NEW_SCATTER: Use new (better) scattering files from Wise.
;   SKIP_SETUP: Do not change marx_commands in the directory,
;               RUN and POST_PROCESS and IMAGE as usual.
;   MARX_SIM_DIR: full path to simulation directory
;
; OTHER SETUP things
; Hardwired location for the simulation results (see in code below.)
; The simulations are stored in this directory under
; the my_obs_id and then an iteration number, e.g., :
;   marx is run in:     marx_sim_dir/H-HAS-EA-8.003/i0/
;   creating output in: marx_sim_dir/H-HAS-EA-8.003/i0/marx_out/
; Note a distribution version of marx.par is assumed to be
; in marx_sim_dir also.

;
if n_elements(my_obs_id) EQ 0 then my_obs_id = 'O-HAS-Capella-3.001'

; read in the obsdb
obsdb = rdb_read('obsdb.rdb')

; find this observation
this_obs = where(STRPOS(obsdb.obs_id,my_obs_id) GE 0, n_found)
if n_found EQ 0 then begin
  print, ''
  print, " Didn't find observation in ObsDB"
  print, ''
  RETURN
end
this_obs = this_obs(0)
obs_id = obsdb(this_obs).obs_id
obs_name = obsdb(this_obs).obs_name
obs_spect = obsdb(this_obs).obs_spect
obs_flux = obsdb(this_obs).obs_flux
obs_time = obsdb(this_obs).obs_time
obs_defocus = obsdb(this_obs).obs_defocus
obs_grating = obsdb(this_obs).obs_grating
obs_detector = obsdb(this_obs).obs_detector

print, ''
print, ' flight_sim for : '+ my_obs_id
print, ''
print, obsdb(this_obs)
print, ''

; OTHER SETUP values, etc.
;- - - - - - - - - - - - - - - - - - - - -
; Directory for output.  marx.par must be here as well.
; Phase 1 simulations:
;;marx_sim_dir = '/nfs/spectra/d8/MARX'
; Phase 2 and Flight simulations:
if n_elements(marx_sim_dir) EQ 0 then marx_sim_dir = '/nfs/spectra/d7/MARX'
; Directory for spectra:
marx_spect_dir = '/nfs/wiwaxia/d4/HETG/GTO_Plan/Marx_spectra'

; MARX distribution directory
marx_dist_dir = '/nfs/wiwaxia/d4/ASC/src/marx-dist-devel'

; "Calibration values":
;
; Flight focal length
flight_foc_len = 10.061 ; m
source_distance = 0  ; m, 0 = infinite

; Add detector and aspect blur: post process uses ray positions
; not FITS files...
; Put ACIS/HRC pixel blur and aspect in here
;  ACIS pixel blur ~ 0.29*24 um = 7.0 um --> 0.14 arc sec rms 1D
;   use 0.14/0.7 = 0.2 arc sec for 2D r rms blur contribution
;  HRC blur is comparable to ACIS...
;  Aspect = 0.34/2.  ; rms radius
source_sigma = SQRT(0.2^2 + 0.17^2)  ; arc seconds rms of 2D gaussian

; Detector and Grating parameters needed for Flight MARX
marx_letg_rs = 8633.92
marx_hetg_rs = 8633.58
;
; Operation paramters
;
numrays_default = 100000 
;- - - - - - - - - - - - - - - - - - - - -

; Make sure the top-level (OBS_ID) directory is present...
found = findfile(marx_sim_dir+'/'+obs_id+'/',COUNT=nfound)
if nfound EQ 0 then SPAWN, 'mkdir '+marx_sim_dir+'/'+obs_id+'/'

; Only one iteration so far...
it_start = 0
it_stop = 0
it_str = STRARR(1)
it_str = ['i0']

if NOT(KEYWORD_SET(SKIP_SETUP)) then begin
; - - - - - inelegant skipping of setup to allow user to run marx then post process ---
; Make the directories and put nominal marx.par in them
for it=it_start,it_stop do begin
  ; iteration directories
  it_dir = marx_sim_dir+'/'+obs_id+'/'+it_str(it)+'/'
  found = findfile(it_dir,COUNT=nfound)
  if nfound EQ 0 then begin
    SPAWN, 'mkdir '+it_dir
    print, '  created dir: ' + it_dir
  end else begin
    print, '   output dir: ' + it_dir
  end
  ; put a copy of marx.par in there
  SPAWN, 'cp '+marx_sim_dir+'/marx.par '+it_dir
end
print, ''


; Put the spectrum file in the directory
for it=it_start,it_stop do begin
  it_dir = marx_sim_dir+'/'+obs_id+'/'+it_str(it)+'/'
  SPAWN, 'cp ' + marx_spect_dir + '/' + obs_spect + ' ' + it_dir
end

; Write a file "marx_commands" in each iteration
; directory.  This file will be sourced to do
; the actual simulation. 
for it=it_start,it_stop do begin
  it_dir = marx_sim_dir+'/'+obs_id+'/'+it_str(it)+'/'
  OPENW, mcunit, it_dir+'marx_commands', /GET
  print, '   creating '+it_dir+'marx_commands'
  printf, mcunit, '# marx_commands for the simulation'
  printf, mcunit, '# created by ~dd/idl/marx/flight_sim.pro ' + SYSTIME()
  printf, mcunit, '# this file is sourced in a directory that already has:'
  printf, mcunit, '#  '+obs_spect
  printf, mcunit, '#    marx.par'
  printf, mcunit, ''
  printf, mcunit, '# where is marx:'
  printf, mcunit, 'set MARX_DIST_DIR = '+marx_dist_dir
  printf, mcunit, ''
  printf, mcunit, '# where is marx data:'
  printf, mcunit, 'setenv MARX_DATA_DIR $MARX_DIST_DIR/marx/data'
  printf, mcunit, ''
  printf, mcunit, '# modify marx.par parameters as needed for the simulation:'
  printf, mcunit, '# - - - - - - - - - - - - - - - - - -'
  printf, mcunit, ''
  printf, mcunit, '# Simulation Setup and Control'
  if n_elements(numrays) EQ 0 then numrays=numrays_default
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par NumRays='+$
	STRCOMPRESS(STRING(numrays),/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par dNumRays=25000'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par ExposureTime=' + $
		STRCOMPRESS(obs_time,/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par OutputDir="marx_out"'
  printf, mcunit, ''
  printf, mcunit, '# Science Instrument Setup and Control'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par GratingType="'+$
			obs_grating+'"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetectorType="'+$
			obs_detector+'"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetIdeal="'+$
			'no'+'"'

  ; Where is the detector?
  ; Defocus
  dX = obs_defocus
  dY = 0.0
  dZ = 0.0

  ; Set marx det offsets appropriately...
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetOffsetX=' + $
		STRCOMPRESS(dX,/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetOffsetY='+$
		STRCOMPRESS(dY,/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetOffsetZ='+$
		STRCOMPRESS(dZ,/REMOVE_ALL)
  printf, mcunit, ''
  printf, mcunit, '# Source Spectral Parameters'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceFlux=' + $
		STRCOMPRESS(obs_flux,/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SpectrumType="FILE"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SpectrumFile="'+ $
			obs_spect+'"'
  printf, mcunit, ''
  printf, mcunit, '# Source Spatial Parameters'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceType="GAUSS"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par S-GaussSigma='+$
			STRCOMPRESS(STRING(source_sigma),/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceDistance='+$
			STRCOMPRESS(STRING(source_distance),/REMOVE_ALL)
  ; Use source elevation and azimuth values here
  src_elev = 0.0
  src_az = 0.0
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceElevation='+$
			STRCOMPRESS(STRING(src_elev),/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceAzimuth='+$
			STRCOMPRESS(STRING(src_az),/REMOVE_ALL)
  printf, mcunit, ''

  printf, mcunit, '# HRMA Setup'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par HRMA_Yaw=0.0'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par HRMA_Pitch=0.0'
  printf, mcunit, ''

  if KEYWORD_SET(NEW_SCATTER) then begin
    printf, mcunit, '# HRMA New Scatter!'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_P1_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_p1_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_H1_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_h1_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_P3_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_p3_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_H3_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_h3_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_P4_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_p4_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_H4_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_h4_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_P6_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_p6_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par WFold_H6_File='+$
			'"/nfs/apocrypha/d3/wise/asc/marx/scat_h6_M"'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par HRMABlur=0.127'
    printf, mcunit, ''
  end

  printf, mcunit, '# Grating Setup and Control'
  rs_value = marx_hetg_rs
  if obs_grating EQ 'LETG' then rs_value = marx_letg_rs
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par RowlandDiameter='+$
			STRCOMPRESS(STRING(rs_value),/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par UseGratingEffFiles="no"'
  printf, mcunit, ''
  printf, mcunit, '# - - - - - - - - - - - - - - - - - -'
  printf, mcunit, ''
  printf, mcunit, '# finally, execute marx:'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/marx @@marx.par'
  printf, mcunit, ''
  printf, mcunit, 'unset MARX_DIST_DIR'
  printf, mcunit, ''
  CLOSE, mcunit
  ; also write the obs_detector and the dX,dY,dZ offset values to a file
  OPENW, mcunit, it_dir+'det_name_xyz.txt'
  printf, mcunit, obs_detector
  printf, mcunit, dX, dY, dZ
  CLOSE, mcunit
  FREE_LUN, mcunit
end
print, ''

; - - - -  end of inelegant skipping of setup...
end

; Now run the simulations!
if KEYWORD_SET(RUN_MARX) then begin
  cd, '', CURRENT=save_idl_dir
  for it=it_start,it_stop do begin
    it_dir = marx_sim_dir+'/'+obs_id+'/'+it_str(it)+'/'

    print, ' Running MARX for iteration '+it_str(it)
    print, ' . . . . . . . . . . . . start MARX . . . . . . . . . . .'
    cd, it_dir
    SPAWN, 'source marx_commands'
    print, ' . . . . . . . . . . . . end MARX . . . . . . . . . . .'
  end
  cd, save_idl_dir
end

; Post-processing things...
if KEYWORD_SET(POST_PROCESS) OR KEYWORD_SET(IMAGE) then begin
  for it=it_start,it_stop do begin
    it_dir = marx_sim_dir+'/'+obs_id+'/'+it_str(it)+'/'
    ; read in the detector name and dX,dY,dZ
    OPENR, det_unit, it_dir+'det_name_xyz.txt',/GET
    det_name = ''
    readf, det_unit, det_name
    readf, det_unit, dX, dY, dZ
    CLOSE, det_unit
    FREE_LUN, det_unit

    print, ' - - - - - - - - - - - - - - '
    print, det_name, dX, dY, dZ
    detX = read_marx_file(it_dir+'marx_out/ypos.dat') - dY
    detY = read_marx_file(it_dir+'marx_out/zpos.dat') - dZ

    window,0
    ; Detector-specific things...
    CASE det_name OF
      'HRC-I': begin
        det_x_range = 70.0
        det_y_range = 70.0
        im_bin = 6.*0.024
        ; HRC PHA goes 0 to 256, scale it to 10 keV
        detE = 10.0*read_marx_file(it_dir+'marx_out/pha.dat')/256.0
        qe_rdb_file = 'N/A'
      end
      'HRC-S': begin
        det_x_range = 150.0
        det_y_range = 15.0
        im_bin = 12.*0.024
        ; HRC PHA goes 0 to 256, scale it to 10 keV
        detE = 10.0*read_marx_file(it_dir+'marx_out/pha.dat')/256.0
        qe_rdb_file = 'N/A'
      end
      'ACIS-S': begin
        det_x_range = 80.0
        det_y_range = 15.0
        im_bin = 6.*0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        ; blur by 100 eV FWHM
        rand_e = detE
        for ie=0L, n_elements(detE)-1 do rand_e(ie) = g_rand(SEED)
        detE = 0.04255*rand_e + detE
        qe_rdb_file = 'N/A'
      end
      'ACIS-I': begin
        det_x_range = 30.0
        det_y_range = 30.0
        im_bin = 4.*0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        ; blur by 100 eV FWHM
        rand_e = detE
        for ie=0L, n_elements(detE)-1 do rand_e(ie) = g_rand(SEED)
        detE = 0.04255*rand_e + detE
        qe_rdb_file = 'N/A'
      end
    ELSE: begin  ;
      print, ' * Not a flight detector!'
      RETURN
    end
    ENDCASE

    select = where(abs(detX) LT det_x_range AND $
		abs(detY) LT det_y_range)


    detX = detX(select)
    detY = detY(select)
    detE = detE(select)

    ; Save these "level 1" products
    ; Create the mL1 structure to return
    ; mL1 has the tags: DETX, DETY, ENERGY
    one_event = {detx:0.0, dety:0.0, energy:0.0}
    mL1 = REPLICATE(one_event, n_elements(detX))
    ; Fill it
    mL1.detx = detX
    mL1.dety = detY
    mL1.energy = detE
    SAVE, mL1, FILENAME=it_dir+'ml1.idlsav'
    ; Return the number of events
    n_post_events = n_elements(mL1)

    plot, detX, detY, PSYM=3, TITLE=obs_id+': Photons hitting '+det_name, $
	XTITLE='Detector X (mm)', YTITLE='Detector Y (mm)', $
	XRANGE=det_x_range*[-1.,1.], YRANGE=det_y_range*[-1.,1.], $
	XSTYLE=1, YSTYLE=1

    if KEYWORD_SET(IMAGE) then begin    
      xye_image, detX, detY, detE, /LOG, $
	GIF=it_dir+'image.gif', ONE=0.5, XBIN=im_bin, YBIN=im_bin

    end
  end
end

RETURN
END