PRO EBRT
; This IRT TRACE procedure does tracing of many grating facets
; of an HETG-like instrument.  Here, the gratings are mounted on
; the Rowland torus and each is oriented to have its normal along the direction
; of the central x-ray path, i.e., as specified in B. et al.
;
; 4/17/91 dd, 4/19/91 dd Added detector on RC
; 5/31/91 dd, add "input" parameters from error budget parameters
;             and report the result expected by the error budget
; 6/2/91 dd, add source angle correction for XRC =/= X, add ACIS tilt.
; 6/7/91 dd, add dp/p for each grating using g_rand function!
; 8/13/92 dd, remove the middle shell for AXAF-I !

; Set plots white on black
@w_on_b

; Include the irt common statements
@irt_common
@xray_common
; Include the ONE_GRATING common
@og_common
; Include the error budget common
@eb_common

; Double precision
PREC = 2

!MTITLE = ' '

; Set flag for those error budget terms that are included in the ray trace:
eb_rt = 0.0 * eb_rt	; reset all to 0
eb_rt(eb_FFS) = 1.0	; Finite-Facet Size error
eb_rt(eb_fD) = 1.0	; Focus and astigmatism error
eb_rt(eb_tD) = 1.0	; Detector Tilts are included in above term
eb_rt(18) = 1.0		; plate curvature error, i.e., non-Rowland use
eb_rt(eb_dppB) = 1.0	; Facet-to-facet dp/p error
eb_rt(eb_u1PH) = 1.0	; HETG plate decenter, perp to grat lines
eb_rt(eb_u2PH) = 1.0	; HETG plate decenter, parallel to grat lines
eb_rt(eb_vPH) = 1.0	; HETG plate defocus
eb_rt(eb_e1PH) = 1.0	; HETG plate tilt, perp to grat lines
eb_rt(eb_e2PH) = 1.0	; HETG plate tilt, parallel to grat lines


; Update the error budget
eb_update

; Initialize ray-trace parameters, many from error budget values
; Operation related parameters:
if n_elements(PRI) eq 0 then PRI = 0
if n_elements(og_PRI) eq 0 then og_PRI = 0
if n_elements(PV) eq 0 then PV = -1
; "Source" parameters:
og_lambda = eb_ps(eb_Lambda)
if n_elements(og_nrays) eq 0 then og_nrays = 100
; Telescope parameters:
og_fl = eb_ps(eb_F)
; Parameters for many gratings
if n_elements(mg_naz) eq 0 then mg_naz = 17
if n_elements(mg_nshells) eq 0 then mg_nshells = 2
; Parameters for the grating facet (in og_common)
og_per = eb_ps(eb_p)
og_order = -1
og_X = eb_ps(eb_XRC)	; Curvature of Rowland torus
ebXhetg = eb_ps(eb_X)	; HETG to HRMA focus distance
og_L = eb_ps(eb_L)

; Detector parameters:
; offset from RC by eb_ps(eb_fD)

; Report parameters:
PRINT, '     Error Budget driven Ray Trace Program '
PRINT, '           Rowland Torus Geometry '
PRINT, ' '
PRINT, 'Parameters:'
PRINT, ' X-ray wavelength          (eb_Lambda) : ', og_lambda, ' A'
PRINT, ' Number of rays per facet     og_nrays : ', og_nrays
PRINT, ' Telescope focal len.           (eb_F) : ', og_fl, ' mm'
PRINT, ' Grating period                 (eb_p) : ', og_per, ' A'
PRINT, ' Diffraction order            og_order : ', og_order
PRINT, ' HETG to HRMA focus  (on axis)  (eb_X) : ', ebXhetg, ' mm'
PRINT, ' HETG Rowland diameter        (eb_XRC) : ', og_X, ' mm'
PRINT, ' Grating size (square)          (eb_L) : ', og_L, ' mm'
PRINT, ' Detector defocus              (eb_fD) : ', eb_ps(eb_fD),'(+ = towards grating)'
PRINT, ' Number of facets (in azimuth)  mg_naz : ', mg_naz
PRINT, ' Number of shells           mg_nshells : ', mg_nshells
IF eb_ps(0) gt 0.5 then $
;mg_thegs = [0.060345, 0.054415, 0.048595]	; MEG and LETG
mg_thegs = [0.060345, 0.048595]	; MEG and LETG, AXAF-I
IF eb_ps(0) lt 0.5 then $
;mg_thegs = [0.042905, 0.037340, 0.031895]	; HEG
mg_thegs = [0.042905,  0.031895]	; HEG AXAF-I
PRINT, ' Shell Ro s (in mm) : ', mg_thegs*ebXhetg
PRINT, ' '

; The rays for each grating facet are calculated and stored in a file
; named by facet number.  Setup for first facet
mg_Nfile=0

; Go through all the facets

FOR mg_M=1, mg_nshells DO BEGIN
source_theg=mg_thegs(mg_M-1)
; source_theg is the off-axis angle for the grating center VIEWED FROM the
; HRMA focus, ie X = ebXhetg away
; Center of entrance pupil, corrected for ebXhetg - XRC mismatch
  og_theg = ATAN((ebXhetg/og_X)*TAN(source_theg))
print, ' og_theg = ', og_theg
FOR mg_N=1, mg_naz DO BEGIN
; Set azimuthal location of grating
  og_phig = (mg_N-1)*(2.0*!PI/mg_naz)

; The beam is a converging set of rays randomly starting in a square
; area of size 2*L by 2*L at the location of the "mirror".
  OBJP, 0., 0., -og_fl, og_nrays, og_lambda	; come to a point focus
;  OBJCR, 0., 0.012, 0., 0., -og_fl, og_nrays, og_lambda	; circle focus
  og_pupxc = og_fl * tan(source_theg) * sin(og_phig)
  og_pupyc = og_fl * tan(source_theg) * cos(og_phig)
; Select random points in entrance pupil
  og_xpup = og_pupxc + 2.0*og_L*(RANDOMU(SEED,og_nrays)-0.5)
  og_ypup = og_pupyc + 2.0*og_L*(RANDOMU(SEED,og_nrays)-0.5)
  og_zpup = FLTARR(og_nrays)  ; all zeros
  PUPIL, og_xpup, og_ypup, og_zpup

; Go to the on-axis grating location
  DISP, 0., 0., (ebXhetg - og_fl)  ; on-axis grating location
  TRAVEL, 3, 0.0
; Put a single grating here with Rowland parameter og_X = eb_XRC
; Period is random due to dp/p between facets
og_per = eb_ps(eb_p)*(1+eb_Ps(eb_dppB)*1.E-4*g_rand())
; Offset (decenter and defocus) the grating plate:
DISP, eb_Ps(eb_u2PH), eb_Ps(eb_u1PH), eb_Ps(eb_vPH)
; Tilt the grating plate, x=parallel to grat lines, y = perp
ROT, eb_Ps(eb_e2PH)/60., 0., 0.
ROT, 0., eb_Ps(eb_e1PH)/60., 0.
;---------------
; Test Cash's torgrath.pro by using it here instead of ONE_GRATING:
; Rotate so dispersion direction is correct
; ROT, 0., 0., -90.
; HETGFULL
; TORGRATM, 1
; ROT, 0., 0., 90.
;---------------
 ONE_GRATING
;---------------
ROT, 0., -eb_Ps(eb_e1PH)/60., 0.
ROT, -eb_Ps(eb_e2PH)/60., 0., 0.
DISP, -eb_Ps(eb_u2PH), -eb_Ps(eb_u1PH), -eb_Ps(eb_vPH)

; Save the rays from this grating in file Nfile.DAT
  mg_Nfile = mg_Nfile+1
  rmat = float(rmat)
  dmat = float(dmat)
  SAVERAYS, STRING(mg_Nfile)

ENDFOR
ENDFOR

;------
; Reconstruct the rays from all gratings
; Get the rays from the first grating
mg_N = 1
GETRAYS, STRING(mg_N)
; add in the rays from all the others
FOR mg_N = 2, mg_Nfile DO ADDRAYS, STRING(mg_N)


; Go to the Gaussian focus
DISP, 0., 0., ( -ebXhetg )
; Move defocus amount and set this as the z value
DISP, 0., 0., eb_ps(eb_fD)
set_nz, eb_ps(eb_fD)

; Set detector on Rowland circle
; Use the grating Rowland circle:
ebRCccd = eb_ps(eb_XRC)
; Use the HETG-focus Rowland circle:
; ebRCccd = ebXhetg
; For a plane detector calculate offset to hit rowland circle...
og_detect = (ebRCccd/2.)*(1.-COS(2.*ASIN(ABS(og_order)*og_lambda/og_per)))
PRINT, ' Detector plane location    og_detect : ', $
	og_detect,'mm (+ = towards grating)'
DISP, 0., 0., ( og_detect )

; x-y plot of intersection with detector plane including tilt effects
ROT, eb_ps(eb_tD)/60., 0., 0.
MOVEZ, 0.
ROT, -eb_ps(eb_tD)/60., 0., 0.

saverays, 'ebrt_last'

; Report what the error budget thinks is going on
eb_rtguess

RETURN
END