;+
;                             FE25_TEMP
;
;  To overlay FeXXV data and synthetic spectra, scaled to peak of w line
;  with scaling of wavelengths as necessary. You must have run the routines
;  GS and GS_CUR first to set up the internal arrays and to define the data
;  time window. eg.
;
;     IDL>  .run gs
;     IDL>  .run gs_cur        (to specify time window)
;     IDL>  .run fe25_temp     (answer prompts)
;
;           The raw spectrum is first displayed and the user should select
;           the resonance line window.  A gaussian is fitted to determine
;           the w-line position and width.  The position determines the
;           wavelength offset and the width is used to determine how much
;           to broaden the synthetic spectrum by.  Since the broadening is
;           by convolution with a gaussian, the spectral fit is not good in
;           the line wings.
;
;           The synthetic spectra plotted are those for the 
;           temperature (+/-1 MK) which is calculated from the ratio of w 
;           to j line intensities.
;
;           Hard copy can be had by running interactively first and then 
;           switching to hardcopy device and rerunning.
;
;           Note that the marked position of the j-line is used to 
;           stretch/compress the scale of the synthetic spectra to match 
;           the observed spectrum, and its intensity is used to calculate the 
;           line ratio.  
;
;
;  CDP  Updated from original t25_plot (now deleted) June 92
;-
 


;
;  NOTE that gs_ds1 & gs_ds2 are the indices in the GSA array of the selected
;  data sets.  The dset_str.dset array maps this to the data file.
;
;
;  allocate spectral array
;
gsp = fltarr(768)
gssz = size(gsa)
if gs_ds1 gt gs_ds2 then begin
   temp = gs_ds1
   gs_ds1 = gs_ds2
   gs_ds2 = temp
endif
if gs_ds1 lt 0 then begin
   print,'Start data set number invalid.'
   goto, fin
endif
if gs_ds2 gt gssz(1)-1 then begin
   print,'End data set number greater than array size.'
   goto, fin
endif

;
;  check for consistent modeid 
;
if (total(index(gs_ds1:gs_ds2).bcs.modeid-index(gs_ds1).bcs.modeid) ne 0) then begin
    print,'Warning: not all MODEIDs in the selected range are the same.'
endif

;
;   and then pick up number of bins per spectrum (omit duff one at end)
;
chan_r = bcs_chan_rng(index,gs_ds1) 
c1s = chan_r(0)
c1e = chan_r(1)-1
c2s = chan_r(2)
c2e = chan_r(3)-1
c3s = chan_r(4)
c3e = chan_r(5)-1
c4s = chan_r(6)
c4e = chan_r(7)-1
;
;  where in the roadmap array are the limits?
;
ds1 = dset_str(gs_ds1).dset
ds2 = dset_str(gs_ds2).dset

;
;  load spectral array
;
if gs_ds1 eq gs_ds2 then begin
   tot_dgi = roadmap(ds1).dgi*0.125
   if gs_decomp eq 0 then begin
      print,'Decompressing data .....'
      gsp = bcs_decomp(gsa(*,gs_ds1))/tot_dgi
   endif else begin
      gsp = gsa(*,gs_ds1)/tot_dgi
   endelse
endif else begin
   tot_dgi = 0.0
   if gs_decomp eq 0 then  print,'Decompressing data .....'
   for i=gs_ds1,gs_ds2 do begin
      if gs_decomp eq 0 then begin
         gsp = gsp + bcs_decomp(gsa(*,i))
      endif else begin
         gsp = gsp + gsa(*,i)
      endelse
      tot_dgi = tot_dgi + roadmap(dset_str(i).dset).dgi*0.125
   endfor

;  correct for total exposure time
;
   gsp = gsp/tot_dgi
endelse
;
;  get wavelength calibrations
;

if n_elements(x_26) lt 100 then get_26,x_26
if n_elements(x_25) lt 100 then get_25,x_25
if n_elements(x_19) lt 100 then get_19,x_19
if n_elements(x_15) lt 100 then get_15,x_15
;
;  subtract observed continuum
;
co = total(gsp(451:470))/20.0
gsp = gsp - co
!x.title='Bins'
!y.title='Counts/sec'
!mtitle='FeXXV spectrum'
!y.range = [0,1.2*max(gsp(255:511))]

if !d.name ne 'PS' then begin
;
;  plot raw data
;
   f = gsp(255:511)
   plot,f,psym=10
   pmax = max(f)
;
;  mark w line position for shift calculation
;
   xyouts,0.25,0.85,'Use 2 cursor clicks to define window on resonance line:',$
                        /normal,charsize=1.3
   print,'Use cursor to define resonance line window'
   nfit = 0
   cursor,s1,qq,3,/data
   if !d.name eq 'TEK' then r = get_kbrd(0)
   oplot,[s1,s1],[0,1000]
   cursor,s2,qq,3,/data
   if !d.name eq 'TEK' then r = get_kbrd(0)
   oplot,[s2,s2],[0,1000]
   if s1 gt s2 then begin
      temp = s1
      s1 = s2
      s2 = temp
   endif
;
;  do gaussian fit
;
   nfit = [fix(s1),fix(s2)]
   x    = findgen(nfit(1)-nfit(0)+1) + nfit(0)
   yfit = findgen(nfit(1)-nfit(0)+1)
   a    = fltarr(6)
   y = f(nfit(0):nfit(1))
   yfit = gaussfit(x,y,a)
;
   pmax = max(yfit)
;
   oplot,x,yfit,psym=-3,thick=2
   !linetype = 2
   oplot,[a(1),a(1)],[0,1000]
   !linetype = 0
endif
;
;  from the width of the fit to the observed line calculate the amount to
;  broaden the synthetic spectrum
;
wids = 0.000447         ; width in lamda of synthetic w line
wido = a(2) * 0.00354   ;   "    "  "    "  observed line
smear = sqrt(wido*wido - wids*wids)  ; assume quadratic addition
bwid = (smear/0.001)*2.0 - 2.0            ; fix for definition of broad width

;
; load synthetic arrays and broaden, note s25_synth should only be loaded the
; first time in
;
get_syn_25,s25,bwid,s25_synth
;
;  wavelength of fitted line
;
wx = x_25(fix(a(1)+0.5))
;
;  replot data
;
!x.title='Wavelength (A)'
plot,x_25,gsp(255:511),psym=10,thick=2,xrange=[1.84,1.89],xstyle=1

;
;  calculate where the w-line is in observed array and synthetic array
;  remember s25(*,0) = synthetic wavelength array
;
xbin = where( abs(x_25-wx) lt 0.0004)
xbin_syn = where(abs(s25(*,0)-1.8505) lt 0.0002)
binw_syn = fix(total(xbin_syn)/n_elements(xbin_syn))
binw = fix(total(xbin)/n_elements(xbin))

;
;  peak of w line
;
yw = max(gsp(256+xbin))
;
;  mark centre to confirm OK
;
!linetype = 2
oplot,[wx,wx],[0,yw]
!linetype = 0
;
;  hence offset between observed and synthetic
;
fe25_off = wx - 1.8505
;
;  so this is where j-line should be
;
obs_j = 1.866 + fe25_off
;
;  adjust for likely scale factor
obs_j = (obs_j - 1.8505)*0.95 + 1.8505
;
;  calculate approx data level around there
;
xbin = where( abs(x_25-obs_j) lt 0.004)
yj = max(gsp(256+xbin))
;
;  and mark it for info
;
if !d.name ne 'PS' then begin
   oplot,[obs_j,obs_j],[yj+pmax/10.,1000]
   oplot,[obs_j,obs_j-0.0004],[yj+pmax/10.,yj+pmax/7.]
   oplot,[obs_j,obs_j+0.0004],[yj+pmax/10.,yj+pmax/7.]
   xyouts,obs_j+0.001,yj+pmax/5.,' j-line',size=2.0
   xyouts,obs_j+0.001,yj+pmax/8.,' here'
   xyouts,obs_j+0.001,yj+pmax/15.,' somewhere'
endif
;
;  now mark actual j-line postion AND peak (for line ratio)
;
if !d.name ne 'PS' then begin
   print,'Now click cursor at PEAK of j line. (NB both coordinates are used): '
   cursor,jx,jy,3,/data
endif
;
;  observed  spectrum is usually stretched/compressed relative to synthetic
;  so work out scale needed, starting at w line
;
xscale = (jx-wx)/0.0155
line_ratio = pmax/jy
print,string('Wavelength scale (w to j) = ',xscale,format='(a,f5.3)')
;
;  approx calibration of ratio v temperature
;
temp = fix((line_ratio+2.2)/0.21)
;
;  error check
;
if temp lt 10 or temp gt 26 then begin
   print,'ERROR: suggested temperature of ',temp,' is outside range,',$
         format='(a,i3,a)'
   repeat begin
      read,'Give required temperature comparison (10-26 MK): ',temp
   endrep until (temp ge 10) and (temp le 26)
endif
;
;  fix up the wavelength scale for plotting of synthetic spectra
;
sx  = s25(*,0)
sx(binw_syn:719) = sx(binw_syn) + (sx(binw_syn:719)-sx(binw_syn))*xscale

;
;  over plot synthetic spectrum of calculated temperature (and +/- 1 MK)
;
;oplot,sx+fe25_off,s25(*,temp)*max(gsp(256:511))/s25(temp,1),thick=1
;oplot,sx+fe25_off,s25(*,temp+1)*max(gsp(256:511))/s25(temp+1,1),psym=3
;oplot,sx+fe25_off,s25(*,temp-1)*max(gsp(256:511))/s25(temp-1,1),psym=3
oplot,sx+fe25_off,s25(*,temp)*pmax/s25(temp,1),thick=1
oplot,sx+fe25_off,s25(*,temp+1)*pmax/s25(temp+1,1),psym=3
oplot,sx+fe25_off,s25(*,temp-1)*pmax/s25(temp-1,1),psym=3

;
;  label
;
xyouts,1.8425,pmax,' Synthetic spectra'
xyouts,1.8425,0.9*pmax,' are for:'
xyouts,1.8425,0.80*pmax,string(temp-1,temp,temp+1,' MK',format='(3i3,a)')


xyouts,0.0,0.0,$
   'Time: '+fmt_tim(index(gs_ds1))+' to '+fmt_tim(index(gs_ds2)),/normal
break_file,infil,disk,dir,fil,ext
xyouts,0.8,0.0,fil+ext,/normal
fin:
!x.title=' '
!y.title=' '
!mtitle=' '
set_xy

end