PRO xrcf_sim, trw_id_in, XYZ_REAL=xyz_real, $
	RUN_MARX=run_marx, ITER_SELECT=iter_select, NUMRAYS=numrays, $
	POST_PROCESS=post_process, N_POST_EVENTS=n_post_events, IMAGE=image, $
	NO_GRATING = no_grating, NEW_SCATTER = new_scatter, $
	SKIP_SETUP=skip_setup, MARX_SIM_DIR=marx_sim_dir, QEIDEAL=qeideal, $
	HXDA_DET_OFFSET=hxda_det_offset
;
; Procedure to do a MARX simulation of an XRCF measurement.
;
; Multiple iterations can occur based on:
;   - POSITION: multiple exposures due to a .loc file (Phase 1)
;   - SHUTTER: shutter scans
;   - ENERGY: source energy scans
;   - note that EE tests are simulated once, the EE iterations
;     are done as a part of the post-processing (pin-hole selection).
;
; INPUTS:
;     trw_id:  'H-HAS-EA-2.001'
;
; KEYWORDS:
;   XYZ_REAL:  specifies the actual X,Y,Z offsets of the detector,
;              in (mm) from the focus.  Not yet implemented - put
;              in as a reminder to do it!
;   RUN_MARX:
;   ITER_SELECT: if present, then only the single iteration
;                is considered.
;   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...
;   NO_GRATING: performs same simulations except grating=NONE and
;               X,Y,Z locs are set to 0,0,0.  The trw_id is modified
;               to be trw_idNG.
;   NEW_SCATTER: no used.
;   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
;   QEIDEAL: do not apply the HXDA detector QE to get more rays
;            from the raytrace
;   HXDA_DET_OFFSET: additional X offset to add in for non-ideal
;                    location of the HXDA detector focal planes
;                    (in mm, + places detector towards HRMA)
; COMMONS and such
; Use the CMDB as the input for the test parameters.
;
; OTHER SETUP things
; Hardwired location for the simulation results (see in code below.)
; The simulations are stored in this directory under
; the trw_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.
;
; 9/23/98 dd Added primey stage angle of 0.0058 radian for
;            HXDA simulations
; 4/07+/99 dd Modify to work with MARX 2.20p8 release.
; 4/8/99 dd Use Effic files as default for HETG and LETG
; 4/8/99 dd Update the QE files used for HXDS detectors
; - - - - - - - - - - - - - - - -

; - - - - - - - - - - - - - - - -
; Use CMDB information...
@cmdb_common
;
if cmdb_using_cal then begin
@cmdb_cal_integers
end else begin
@cmdb_trw_integers
end

; OTHER SETUP and XRCF CALIBRATION values, etc.
;- - - - - - - - - - - - - - - - - - - - -
; Directory for output.  marx.par must be here as well.
marx_sim_dir = !DDPHASE1SIM

; MARX distribution directory
marx_dist_dir = !DDMARX220

; Locations files:
loc_file_dir = !DDHETGCAL+'/data/Loc_files'

; "Calibration values":
;
; XRCF focal length
xrcf_foc_len = 10.2586 ; m
; Flight focal length
flight_foc_len = 10.0655 ; m
; HRMA Blur for XRCF will be in the marx.par file
; as there are many blur parameters!
; XRCF source distance and nominal spot size for MARX
source_distance = 527.522 + flight_foc_len  ; m
; Will use a DISK model, below.
source_spot = 0.4 ; mm diameter
source_radius = 0.5*source_spot/(1000.0*source_distance) ; radians
source_radius = 3600.0*source_radius/!DTOR
; Rowland Spacing and Detector Offset for MARX 2.20p8:
xrcf_det_x_offset = -194.832 ; mm
marx_letg_rs = 8788.04 + xrcf_det_x_offset
marx_hetg_rs = 8782.80 + xrcf_det_x_offset
; primey rotation angle (about Z axis):
primey_rot = 0.0058  ; radians, from Brad W.
;
; Operation paramters
;
numrays_default = 100000 
;- - - - - - - - - - - - - - - - - - - - -

; Catch a null trw_id:
if n_elements(trw_id_in) NE 1 then begin
  print, ''
  print, " Usage: xrcf_sim, 'H-HAS-EA-2.001' [, /RUN_MARX, ITER=2, " + $
		"NUMRAYS=100000, /POST_PROCESS, N_POST_EVENTS, /IMAGE, /NO_GRATING, /SKIP_SETUP, /QEIDEAL ]"
  print, ''
  RETURN
end

; Find the measurement in the CMDB
this_row = where(cmdb_fields(*,c_id) EQ trw_id_in, nwhere)
if nwhere EQ 1 then begin
  this_row = this_row(0)
end else begin
  print, ' Found '+trw_id_in+STRING(nwhere,FORMAT='(I2)')+' times in CMDB ?!'
  print, ''
  RETURN
end
; Show the CMDB summary info
print, ''
print, '  CMDB summary information:'
cmdb_mini_list, this_row
print, ''

; Set trw_id
trw_id = trw_id_in
if KEYWORD_SET(NO_GRATING) then begin
  trw_id = trw_id_in + 'NG'
  print, '* Simulation with NO_GRATING'
  print, ''
end

; Flag any items not simulated...
; FAM pitch and yaw not simulated
fam_pitch_str = cmdb_fields(this_row,c_five_axis_pitch)
fam_yaw_str = cmdb_fields(this_row,c_five_axis_yaw)
if (fam_pitch_str NE '0' AND fam_pitch_str NE 'N/A') OR $
	(fam_yaw_str NE '0' AND fam_yaw_str NE 'N/A') then begin
  print, '* FAM pitch,yaw ignored: '+fam_pitch_str+', '+fam_yaw_str
  print, '  However, for non-zero HRMA Pitch and Yaw the detector'
  print, '  is pitched, yawed to maintain flight-like simulation.'
  print, ''
end
; FAM iterations are not currently implemented, but the FAM choice
; is printed if non standard.
; What is the five_axis up to?
fam_choice = cmdb_fields(this_row,c_five_axis_choices)
if fam_choice NE 'N/A' AND fam_choice NE 'FIXED' then begin
  print, '* FAM choice is set for: '+fam_choice+'; not simulated!'
  print, ''
end

; Check for iterations...
; Iterations are exclusive so set the type and number of iterations:
; it_type = [ 'NONE' | 'ENERGY' | 'SHUTTER' | 'POSITION' ]
it_type = 'NONE'
it_num = 1
it_str = ['i0']
it_param = 0

; ENERGY
mono_energies = cmdb_mono_energies(this_row)
if n_elements(mono_energies) GT 1 then begin
  it_type = 'ENERGY'
  it_num = n_elements(mono_energies)
  it_str = STRARR(it_num)
  for it=0,it_num-1 do begin
    it_str(it) = 'i'+STRCOMPRESS(STRING(it),/REMOVE_ALL)
  end
  it_param = mono_energies
end
; SHUTTER
; (For a list of all shutter values actually used at XRCF
;  see ~dd/xrcf/shutters_used.txt)
; Shutter iterations are signalled by a SCAN in the quadrant designation:
; e.g., {1,SCAN}, {ALL,SCAN}, {HEG,SCAN}.
; Get the shutter string for this trw_id
shut_str = cmdb_fields(this_row, c_shutters)
; Done if shutters are closed
if STRPOS(shut_str,'CLOSED') GE 0 then begin
  print, ' Shutters are closed - no simulation!'
  print, ''
  RETURN
end
; OK, check for SCAN and, if so, fill the it_param
; with DISCRETE shutter commands as appropriate.
if STRPOS(shut_str,'SCAN') GE 0 then begin
  it_type = 'SHUTTER'
  it_num = 4  ; SCAN is always the four quadrants
  it_str = STRARR(it_num)
  for it=0,it_num-1 do begin
    it_str(it) = 'i'+STRCOMPRESS(STRING(it),/REMOVE_ALL)
  end
  it_param = STRARR(4)
  ; Fill it_param with the DISCRETE command
  ; The SCAN is done in the order T,N,B,S
  ; Determine the mirror set scanned:
  ; find the ","
  comma_pos = STRPOS(shut_str,',')
  ; start at 1 to remove the "{"
  shell_str = STRMID(shut_str, 1, comma_pos-1)
  CASE shell_str OF
    '1': mir_str = '1'
    '3': mir_str = '3'
    '4': mir_str = '4'
    '6': mir_str = '6'
    'ALL': mir_str = '1346'
    'HEG': mir_str = '46'
    'MEG': mir_str = '13'
    ELSE: begin
      print, ' Unrecognized shell part of '+shut_str
      print, ''
      RETURN
    end
  ENDCASE
  ; FIll them in
  it_param(0) = '{DISCRETE,T'+mir_str+',N-,B-,S-}'
  it_param(1) = '{DISCRETE,T-,N'+mir_str+',B-,S-}'
  it_param(2) = '{DISCRETE,T-,N-,B'+mir_str+',S-}'
  it_param(3) = '{DISCRETE,T-,N-,B-,S'+mir_str+'}'
end
; POSITION
; Iterations over detector positions are determined
; with a locations file, c_loclist_file:
locfile = cmdb_fields(this_row, c_loclist_file)
if locfile NE 'N/A' then begin
  ; Read-in the Loc file
  xrcf_loc_file_read, trw_id_in, Xloc, Yloc, Zloc, $
	LOC_DIR=loc_file_dir
  ; if more than one point then it's an iteration
  ; (note that single point loc files were used
  ;  quite often by gratings folks!)
  if n_elements(Xloc) GT 1 then begin
    it_type = 'POSITION'
    it_num =  n_elements(Xloc)
    it_str = STRARR(it_num)
    for it=0,it_num-1 do begin
      it_str(it) = 'i'+STRCOMPRESS(STRING(it),/REMOVE_ALL)
    end
    ; Fill it_param with X,Y,Z locations w.r.t. best focus
    ; comma delimited:
    it_param = STRARR(it_num)
    xyzfmt = '(F10.3)'
    for it=0,it_num-1 do begin
      it_param(it) = STRING(Xloc(it),FORMAT=xyzfmt) + ',' + $
			STRING(Yloc(it),FORMAT=xyzfmt) + ',' + $
			STRING(Zloc(it),FORMAT=xyzfmt)
    end
  end
end

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

; Write a trw_id_iters.rdb file:
if it_type EQ 'ENERGY' then begin
  one_it = {type:'', str:'', param:0.0}
end else begin
  one_it = {type:'', str:'', param:''}
end
its_rdb = REPLICATE(one_it, it_num)
for in=0,it_num-1 do begin
  its_rdb(in).type = it_type
  its_rdb(in).str = it_str(in)
  its_rdb(in).param = it_param(in)
end
rdb_write, marx_sim_dir+'/'+trw_id+'/'+trw_id+'_iters.rdb', its_rdb
print, ''

; Summarize the iterations, if any:
if it_num GT 1 then begin
  print,' Simulation has multiple iterations in '+it_type
  for it=0, it_num-1 do begin
    print, '     ' + it_str(it)+' : ', it_param(it)
  end
end else begin
  print,' Simulation has a single iteration, i0'
end


; The above setup activities occur for all iterations
; The directory setup, marx.par creation, running of MARX
; and any post processing can be
; done on a single iteration as selected by ITER_SELECT
; keyword.
if N_ELEMENTS(ITER_SELECT) GT 0 then begin
  it_start = ( (iter_select > 0) < (it_num-1) )
  it_stop = ( (iter_select > 0) < (it_num-1) )
end else begin
  it_start = 0
  it_stop = it_num-1
end

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+'/'+trw_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, ''

; Create the spectrum file(s)
; They have the name trw_id.spec
; If the it_start one exists then skip this...
found = findfile(marx_sim_dir+'/'+trw_id+'/'+ $
		it_str(it_start)+'/'+trw_id_in+'.spec',COUNT=nfound)
; If spectra have changed, etc., it can be useful to force a new
; spectral calculation:
;;nfound = 0  ; force new spectrum calc...
if nfound EQ 0 then begin
  if it_type EQ 'ENERGY' then begin
    ; Run xss_sim for each selected iteration
    cd, '', CURRENT=save_idl_dir
    for it=it_start,it_stop do begin
      ; iteration directories
      it_dir = marx_sim_dir+'/'+trw_id+'/'+it_str(it)+'/'
      cd, it_dir
      xss_sim, this_row+1, ITER = it
    end
    cd, save_idl_dir
  end else begin
    ; Run xss_sim for it_start and copy to rest
    cd, '', CURRENT=save_idl_dir
    it = it_start
    it_dir = marx_sim_dir+'/'+trw_id+'/'+it_str(it)+'/'
    cd, it_dir
    xss_sim, this_row+1, ITER = it
    it0_dir = it_dir
    cd, save_idl_dir
    for it=it_start+1,it_stop do begin
      it_dir = marx_sim_dir+'/'+trw_id+'/'+it_str(it)+'/'
      SPAWN, 'cp '+it0_dir + trw_id+'.* '+it_dir
    end
  end
end ; making spectrum files if needed

; 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+'/'+trw_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/xrcf_sim.pro ' + SYSTIME()
  printf, mcunit, '# this file is sourced in a directory that already has:'
  printf, mcunit, '#  '+trw_id+'.spec'
  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=0.0'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par OutputDir="marx_out"'
  printf, mcunit, ''
  printf, mcunit, '# Science Instrument Setup and Control'
  grat_string = cmdb_fields(this_row,c_grating)
  grat_string = STRCOMPRESS(grat_string,/REMOVE_ALL)
  if KEYWORD_SET(NO_GRATING) then grat_string = 'NONE'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par GratingType="'+$
			grat_string+'"'
  cmdb_det = cmdb_fields(this_row,c_detector)
  CASE cmdb_det OF
    'FPC2': begin
              marx_det_name = 'HRC-I'
              det_name = 'FPC'
              det_ideal = 'yes'
            end
    'HSI': begin
              marx_det_name = 'HRC-I'
              det_name = 'HSI'
              det_ideal = 'yes'
            end
    'SSD': begin
              marx_det_name = 'HRC-I'
              det_name = 'SSD'
              det_ideal = 'yes'
            end
    'ACIS,2C0': begin
              marx_det_name = 'HRC-I'
              det_name = '2C0'
              det_ideal = 'yes'
            end
    'ACIS,2C1': begin
              marx_det_name = 'HRC-I'
              det_name = '2C1'
              det_ideal = 'yes'
            end
    'HRC,I': begin
              marx_det_name = 'HRC-I'
              det_name = 'HRC-I'
              det_ideal = 'no'
            end
    'HRC,S': begin
              marx_det_name = 'HRC-S'
              det_name = 'HRC-S'
              det_ideal = 'no'
            end
    ELSE: begin
      if STRPOS(cmdb_det,'ACIS,I') GE 0 then begin
        marx_det_name = 'ACIS-I'
        det_name = 'ACIS-I'
        det_ideal = 'no'
      end 
      if STRPOS(cmdb_det,'ACIS,S') GE 0 then begin
        marx_det_name = 'ACIS-S'
        det_name = 'ACIS-S'
        det_ideal = 'no'
      end 
    end
  ENDCASE
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetectorType="'+$
			marx_det_name+'"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetIdeal="'+$
			det_ideal+'"'
  ; Where is the detector?
  if locfile NE 'N/A' then begin
    ; already read in the loc file...
    if it_type EQ 'POSITION' then begin
      dX = Xloc(it)
      dY = Yloc(it)
      dZ = Zloc(it)
    end else begin
      dX = Xloc(0)
      dY = Yloc(0)
      dZ = Zloc(0)
    end
  end else begin
    ; get the location from CMDB info...
    ; Y, Z, and defocus
    if cmdb_using_cal then begin
      ystr=cmdb_fields(this_row,c_offset_y)
      zstr=cmdb_fields(this_row,c_offset_z)
      xstr=cmdb_fields(this_row,c_defocus)
    end else begin
      if STRPOS(cmdb_fields(this_row,c_measurement_config),'HXDA') $
							GE 0 then begin
        ; Using HXDA
        ystr=cmdb_fields(this_row,c_HXDA_Y)
        zstr=cmdb_fields(this_row,c_HXDA_Z)
        xstr=cmdb_fields(this_row,c_HXDA_defocus)
      end else begin
        ; Using FAM etc.
        ystr=cmdb_fields(this_row,c_five_axis_y)
        zstr=cmdb_fields(this_row,c_five_axis_Z)
        xstr=cmdb_fields(this_row,c_five_axis_X)
      end
    end
    READS, ystr, dY
    READS, zstr, dZ
    READS, xstr, dX
  end
  ; dX, dY, and dZ have mm w.r.t. on-axis best focus in XRCF coord.s.

  ; But wait! If the HXDA was used then we need to include (simulate)
  ; the primey rotation effect on the X position of the detector:
  if STRPOS(cmdb_fields(this_row,c_measurement_config),'HXDA') $
	GE 0 then begin
    if n_elements(HXDA_DET_OFFSET) EQ 0 then HXDA_DET_OFFSET = 0.0
    dX = dX - primey_rot * dY + HXDA_DET_OFFSET
  end

  ; For the no-grating case these should all be 0:
  if KEYWORD_SET(NO_GRATING) then begin
    dX = 0.*dX
    dY = 0.*dY
    dZ = 0.*dZ
  end

  ; Convert detector coord.s to MARX coord.s
  dX = dX
  dY = -1.0 * dY
  dZ = -1.0 * dZ
  ; Now detector offset coord.s are in MARX/flight system.

  ; Check and adjust for non-zero HRMA pitch or yaw
  ; This involves i) setting the source off-axis
  ;              ii) modifying the detector offsets
  ; Get the HRMA off-axis info
  hrma_pitch_str = cmdb_fields(this_row,c_HRMA_pitch)
  hrma_yaw_str = cmdb_fields(this_row,c_HRMA_yaw)
  ; From XRCF and MARX geometry we get that:
  ;    Source elevation should be MINUS one times the HRMA pitch
  ;    Source azimuth should be PLUS one times the HRMA yaw
  src_elev = 0.0
  src_az = 0.0
  if hrma_pitch_str NE '0' OR hrma_yaw_str NE '0' then begin
    ; i) Set the source elevation and azimuth.
    READS, hrma_pitch_str, src_elev
    src_elev = -1.*src_elev
    READS, hrma_yaw_str, src_az
    print, '     HRMA pitch, yaw --> Source elevation, azimuth: ', $
		src_elev, src_az
    ; ii) Modify the detector offsets...
    ; Use HRMA xrcf focal length, xrcf_foc_len (meters), to produce
    ; the equivalent displacement the detector has moved to track
    ; the ridgid HRMA-Det motion
    marx_dY = xrcf_foc_len*1000.0 * !DTOR*src_az/60.0
    marx_dZ = xrcf_foc_len*1000.0 * !DTOR*src_elev/60.0
    ; The difference is additional motion needed to have
    ; marx agree with CMDB/XRCF specification:
    dY = dY - marx_dY     
    dZ = dZ - marx_dZ     
    print, '     Det location (MARX coords):', dX, dY, dZ
  end

  ; Set marx det offsets appropriately...
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par DetOffsetX=' + $
		STRCOMPRESS(dX+xrcf_det_x_offset,/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=0.0'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SpectrumType="FILE"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SpectrumFile="'+ $
			trw_id_in+'.spec"'
  printf, mcunit, ''
  printf, mcunit, '# Source Spatial Parameters'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceType="DISK"'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par S-DiskTheta0='+$
			'0.0'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par S-DiskTheta1='+$
			STRCOMPRESS(STRING(source_radius),/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
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceOffsetZ='+$
			STRCOMPRESS(STRING(src_elev),/REMOVE_ALL)
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par SourceOffsetY='+$
			STRCOMPRESS(STRING(src_az),/REMOVE_ALL)
  printf, mcunit, ''
  printf, mcunit, '# XRCF Shutter Control'
  ; SCANs have been converted to DISCRETEs so no need to check for
  ; SCAN.  CLOSED has also been checked, so that leaves DISCRETES and
  ; the following: {1,ALL},{3,ALL},{4,ALL},{6,ALL},{ALL,ALL},
  ;   {HEG,ALL}, {MEG,ALL}
  ; 
  ; If this is a SHUTTER iteration get the string from it_param, 
  ; otherwise from the shutter string:
  if it_type EQ 'SHUTTER' then begin
    shut_str = it_param(it)
  end else begin
    shut_str = cmdb_fields(this_row, c_shutters)
  end
  ; Convert to DISCRETEs:
  CASE shut_str OF
    '{1,ALL}': shut_str = '{DISCRETE,T1,N1,B1,S1}'
    '{3,ALL}': shut_str = '{DISCRETE,T3,N3,B3,S3}'
    '{4,ALL}': shut_str = '{DISCRETE,T4,N4,B4,S4}'
    '{6,ALL}': shut_str = '{DISCRETE,T6,N6,B6,S6}'
    '{ALL,ALL}': shut_str = '{DISCRETE,T1346,N1346,B1346,S1346}'
    '{HEG,ALL}': shut_str = '{DISCRETE,T46,N46,B46,S46}'
    '{MEG,ALL}': shut_str = '{DISCRETE,T13,N13,B13,S13}'
    ELSE:
  ENDCASE
  ; Now convert DISCRETE to MARX code...
  ; Because the XRCF coord.s are upside down w.r.t. flight(MARX)
  ; coord.s, the shutters are "rotated" and mapped as:
  ;  updated 8/24/98 dd
  ;  ShuttersX="ABCD" where:
  ;   A = 'Bottom' (+Z MARX, -Z XRCF)
  ;   B = 'North'  (+Y MARX, -Y XRCF)
  ;   C = 'Top'    (-Z MARX, +Z XRCF)
  ;   D = 'South'  (-Y MARX, +Y XRCF)

  ; Make sure it is a DISCRETE...
  if STRPOS(shut_str, '{DISCRETE') LT 0 then begin
    print, ' * Unknown shutter string: '+shut_str
    print, ' --> All shutters left open'
  end else begin
    ; Do the conversion...shell by shell
    pieces = STR_SEP(shut_str,',')
    tnbs = pieces(1:4)
    ; These next two are not correct!
    ; neworder = [2,3,0,1]
    ; bstn = tnbs(neworder)
    ; correct versions:
    neworder = [2,1,0,3]
    bnts = tnbs(neworder)
    oneifmt = '(I1)'
    ; now query each shell's shutters...
    ; (could make this a loop over shells)
    shell = '1'
    shut_code = 1*(STRPOS(bnts,shell) EQ -1)    ; 1* is to convert to integer
    shut_code_str = STRING(shut_code(0),FORMAT=oneifmt)
    for is=1,3 do shut_code_str = shut_code_str + $
				STRING(shut_code(is),FORMAT=oneifmt)
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par Shutters1="'+$
				shut_code_str+'"'
    shell = '3'
    shut_code = 1*(STRPOS(bnts,shell) EQ -1)    ; 1* is to convert to integer
    shut_code_str = STRING(shut_code(0),FORMAT=oneifmt)
    for is=1,3 do shut_code_str = shut_code_str + $
				STRING(shut_code(is),FORMAT=oneifmt)
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par Shutters3="'+$
				shut_code_str+'"'
    shell = '4'
    shut_code = 1*(STRPOS(bnts,shell) EQ -1)    ; 1* is to convert to integer
    shut_code_str = STRING(shut_code(0),FORMAT=oneifmt)
    for is=1,3 do shut_code_str = shut_code_str + $
				STRING(shut_code(is),FORMAT=oneifmt)
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par Shutters4="'+$
				shut_code_str+'"'
    shell = '6'
    shut_code = 1*(STRPOS(bnts,shell) EQ -1)    ; 1* is to convert to integer
    shut_code_str = STRING(shut_code(0),FORMAT=oneifmt)
    for is=1,3 do shut_code_str = shut_code_str + $
				STRING(shut_code(is),FORMAT=oneifmt)
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par Shutters6="'+$
				shut_code_str+'"'
    printf, mcunit, ''
  end
  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, ''

  printf, mcunit, '# Grating Setup and Control'
  rs_value = marx_hetg_rs
  if grat_string 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="yes"'

  printf, mcunit, ''
  ; This is set to zero when HRC-I is being used to model an HXDS instrument
  if STRPOS(cmdb_det,'HRC') LT 0 then begin
    printf, mcunit, '# HRC Model Parameters'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/pset marx.par HRC-I-BlurSigma'+$
				'=0.0'
    printf, mcunit, ''
  end
  printf, mcunit, '# - - - - - - - - - - - - - - - - - -'
  printf, mcunit, ''
  printf, mcunit, '# finally, execute marx:'
  printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/marx @@marx.par'
  printf, mcunit, ''
  ; For ACIS and HRC do marx2fits also
  if STRPOS(det_name, '-') GT 0 then begin
    printf, mcunit, ''
    printf, mcunit, '# Do marx2fits for this flight detector:'
    printf, mcunit, '$MARX_DIST_DIR/bin/$ARCH/marx2fits marx_out '+ $
	'marx_out.fits'
    printf, mcunit, ''
  end
  printf, mcunit, 'unset MARX_DIST_DIR'
  printf, mcunit, ''
  CLOSE, mcunit
  ; also write the det_name and the dX,dY,dZ offset values to a file
  OPENW, mcunit, it_dir+'det_name_xyz.txt'
  printf, mcunit, det_name
  printf, mcunit, dX, dY, dZ
  CLOSE, mcunit
  FREE_LUN, mcunit
end
print, ''

; - - - -  end of inelegant....
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+'/'+trw_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+'/'+trw_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...
    hxds_cip_dir = !DDHXDSDIR + '/'
    CASE det_name OF
      'FPC': begin
        det_x_range = 6.0
        det_y_range = 6.0
        facility_coords = 1
        im_bin = 0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        qe_rdb_file = 'fpc_x2.fitQe_N19980902.rdb'
        det_qe = rdb_read(hxds_cip_dir+qe_rdb_file)
        ; energy is in keV
      end
      'SSD': begin
        det_x_range = 6.0
        det_y_range = 6.0
        facility_coords = 1
        im_bin = 0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        qe_rdb_file = 'ssd_x_qe_N19990208.rdb'
        det_qe = rdb_read(hxds_cip_dir+qe_rdb_file)
        ; Energy is in eV, convert to keV
        det_qe.energy = det_qe.energy/1000.0
      end
      'HSI': begin
        det_x_range = 10.0
        det_y_range = 10.0
        facility_coords = 1
        im_bin = 0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        qe_rdb_file = 'hsi_qe_N19961106.rdb'
        det_qe = rdb_read(hxds_cip_dir+qe_rdb_file)
      end
      'HRC-I': begin
        det_x_range = 70.0
        det_y_range = 70.0
        facility_coords = 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
        facility_coords = 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
        facility_coords = 0
        im_bin = 6.*0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        qe_rdb_file = 'N/A'
      end
      'ACIS-I': begin
        det_x_range = 30.0
        det_y_range = 30.0
        facility_coords = 0
        im_bin = 4.*0.024
        detE = read_marx_file(it_dir+'marx_out/energy.dat')
        qe_rdb_file = 'N/A'
      end
    ELSE: begin  ; 2C
      det_x_range = 10.0
      det_y_range = 10.0
      facility_coords = 1
      im_bin = 0.024
      detE = read_marx_file(it_dir+'marx_out/energy.dat')
      qe_rdb_file = 'N/A'  ; will need file for the 2 2C's
    end
    ENDCASE

    ; Window out the '[H|L]AS-MC-3.' zero order +/- 1.5 mm:
    if STRPOS(trw_id,'HAS-MC-3.') GE 0 OR $
	STRPOS(trw_id,'LAS-MC-3.') GE 0 then begin
      select = where( abs(detX) LT det_x_range AND $
		abs(detY) LT det_y_range AND $
		abs(detX) GT 1.5)
    end else begin
      select = where(abs(detX) LT det_x_range AND $
		abs(detY) LT det_y_range)
    end
    ; and 'MC-20.001
    if STRPOS(trw_id,'HAS-MC-20.001') GE 0 then begin
      print, ' * MC-20.001 : windowing out zero-order...'
      select = where( (abs(detX) LT det_x_range) AND $
		(abs(detY) LT det_y_range) AND $
		(abs(detX) GT 1.2 OR abs(detY-1.0) GT 1.0) )
    end

    if facility_coords then begin
      detX = -1.*detX(select)
      detY = -1.*detY(select)
      detE = detE(select)
    end else begin
      detX = detX(select)
      detY = detY(select)
      detE = detE(select)
    end

    ; Apply detector QE file if required (HXDS detectors)
    if qe_rdb_file NE 'N/A' AND NOT(KEYWORD_SET(QEIDEAL))then begin
      ; det_qe has tags energy and qe
      ; assign a uniform random variable to each photon
      urv = RANDOMU(SEED, n_elements(detE))
      ; calc the QE at the energy of each photon
      qe_for_phot = INTERPOL(det_qe.qe, det_qe.energy, detE)
      ; select only those photons where QE GT random variable
      select = where(qe_for_phot GT urv)
      detX = detX(select)
      detY = detY(select)
      detE = detE(select)
    end

    ; 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)

    if facility_coords then begin
      plot, detX, detY, PSYM=3, TITLE=trw_id+': Photons hitting '+det_name, $
	XTITLE='Y_XRCF (mm)', YTITLE='Z_XRCF (mm)', $
	XRANGE=det_x_range*[-1.,1.], YRANGE=det_y_range*[-1.,1.], $
	XSTYLE=1, YSTYLE=1
    end else begin
      plot, detX, detY, PSYM=3, TITLE=trw_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
    end

    if det_name EQ 'HSI' then begin
      hsi_mask, x_mask, y_mask
      oplot, x_mask, y_mask, LINESTYLE=1
    end

    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

      ; Put the XRCF data and the MARX data in the same image
      if det_name EQ 'ACIS-S' then begin
        marx_compare, trw_id, /REST, /IMAGE, GIF=it_dir+trw_id+'_compare.gif'
      end
      if det_name EQ 'HRC-I' then begin
        marx_compare, trw_id, /REST, /IMAGE, GIF=it_dir+trw_id+'_compare.gif'
      end
    end
  end
end

RETURN
END