;+
;
; Project   : s3drs reconstruction software
;                   
; Name      : (various test routines)
;               
; Purpose   : part of the s3drs test kit
;               
; Explanation: (archived here for completeness, not intended for users use)
;               
; Use       : (self-documented, see notes in code)
;    
; Written     : Sandy Antunes, NRC, 2005-2009
;               
;-            
; test script that takes a generic SECCHI file,
; then inserts data representing different test
; cases of '1DN' and '100DN' to compare the results
; of secchi_prep and extract numerics, in order to
; better help us figure out photons and noise.


; Since secchi_prep returns DN, we then create our fake 0.5 noise
; files using a 'clean' header with bscale=1, bzero=0 in order to
; reliably create something that will pass secchi_prep processing.


 on_error,2

  plot=0; 1 = yes, 0 = no

; ###########################
; # PREPARATION TIME
; # ingest a SECCHI file, to use its header for all future work
 fname='cor105.fts'
 secchi_prep,fname,/write_fts; just to make sure it works
; # this is the default output file name for this input file
 outfname='20061106_131400_15c1A.fts'
 if (file_exist(outfname) eq 0) then print,'wrong output file name ',outfname

; # now run 'secchi_prep' on it, saving the original as well
 fname1DNa='1DNcor1_L05.fts'
 fname1DNb='1DNcor1_L10.fts'

; hdr=fitshead2struct(headfits(fname))
; data=readfits(fname)
 data_orig=sccreadfits(fname,hdr)
 help,data_orig

 ; soho offsetbias = 332; LASCO C1=332, C2=470, C3=319, 
 biasmean = hdr.biasmean

 ; these two should be 0, 1 in any real SECCHI header
 hdr.bzero=0
 hdr.bscale=1
 bzero = hdr.bzero
 bscale = hdr.bscale

; ############################
; # Make our fake data sets for 1DN and 100DN, to run under
; # secchi_prep to see what secchi_prep does to different levels
; # of incident signal.
; # FIRST make and save 1DN file

;; data = fix((data * 0) + 1 + biasmean)
 rawdata1DN = uint( ((data_orig * 0) + 1 + biasmean - bzero)/bscale )
 help,rawdata1DN
 hdranon=create_struct(hdr,name='')
 hdranon=rep_tag_name(hdranon,'EXTEND','SIMPLE')
 hdranon=rep_tag_value(hdranon,'1','SIMPLE')
 hdrstr = struct2fitshdr(hdranon)
 writefits,fname1DNa,rawdata1DN,hdrstr
; # Now prep it
 secchi_prep,fname1DNa,/write_fts
 rename,outfname,fname1DNb

; # AND also make a 128DN file

 fname128DNa='128DNcor1_L05.fts'
 fname128DNb='128DNcor1_L10.fts'
; data = fix((data * 0) + 128 + biasmean)
 rawdata128DN = uint( ((data_orig * 0) + 128 + biasmean - bzero)/bscale )
; writefits,fname128DNa,rawdata128DN,struct2fitshdr(hdr)
 writefits,fname128DNa,rawdata128DN,hdrstr
 secchi_prep,fname128DNa,/write_fts
 rename,outfname,fname128DNb

; # so files are now fname1DNa/b:   1DNcor1_L05.fts, 1DNcor1_L10.fts
; #                 fname128DNa/b: 128DNcor1_L05.fts, 128DNcor1_L10.fts
; # also original flat field file cor1flatfield.fts
  flatfield='cor1flatfield.fts'

; # now read in the different data sets
  data_1DN_05=sccreadfits(fname1DNa)
  data_1DN_10=sccreadfits(fname1DNb)
  data_128DN_05=sccreadfits(fname128DNa)
  data_128DN_10=sccreadfits(fname128DNb)
  data_flattemp=readfits(flatfield); or should I use sccreadfits?
;  data_flat=data_flattemp[102:*,102:*]
;  data_flat=data_flattemp[50:2097,50:2097]
  data_flat=data_flattemp[hdr.P1COL:hdr.P2COL:,hdr.P1ROW:hdr.P2ROW]
; # make smaller 128x128 versions of it
  data_1DN_05c=congrid(data_1DN_05,128,128)
  data_1DN_10c=congrid(data_1DN_10,128,128)
  data_128DN_05c=congrid(data_128DN_05,128,128)
  data_128DN_10c=congrid(data_128DN_10,128,128)
  data_flatc=congrid(data_flat,128,128)

; ############################
; # now make an artificial file equivalent to 1 DN
; # for both 1 sec exptime and 100 sec exptime, so we can
; # see how secchi_prep differs for the same data accumulated
; # over different exposure times
  dn2msb=7.1E-11
;;  data_1DN_art=((data*0)+1+biasmean)
;  data_1DN_art = uint( ((data_orig * 0) + 1 + biasmean - bzero)/bscale )
;  hdr.exptime=1
;  fname1DNa1sec='1sec05.fts'
;  writefits,fname1DNa1sec,data_1DN_art,struct2fitshdr(hdr)
;  writefits,fname1DNa1sec,rawdata1DN,hdrstr

 hdranonB=rep_tag_value(hdranon,100,'EXPTIME')
 hdrstrB = struct2fitshdr(hdranonB)
  fname1DN100secL05='100sec05.fts'
  fname1DN100secL10='100sec10.fts'
;  writefits,fname1DNa100sec,data_1DN_art,struct2fitshdr(hdr)
  writefits,fname1DN100secL05,rawdata1DN,hdrstrB

  ; # prep these fake SECCHI/SOHO data up to level 1
;  secchi_prep,fname1DNa1sec
;  rename,outfname,fname1DN1secL10
  secchi_prep,fname1DN100secL05,/write_fts
  rename,outfname,fname1DN100secL10

; # and create a dataset of what they should end up as, i.e. fake secchi_prep
  data_msb= ((data_orig*0)+1) * dn2msb * data_flat

; ###########################
;# Now make a fake SECCHI file with inserted LASCO data, and prep it to Level 1

;  sohoin='/home/antunes/data_old/22003802.fts'
;  sohoout05='hybrid05.fts'
;  sohoout10='hybrid10.fts'
;  sohodata=congrid(readfits(sohoin),2048,2048)
;  writefits,sohoout05,sohodata,hdrstr
;  secchi_prep,sohoout05
;  rename,outfname,sohoout10

print,'biasmean = ',biasmean
print,'1DN 1 sec L0.5 (DN) mean = ',mean(data_1DN_05),$
  ' -biasmean * dn2msb = ',(mean(data_1DN_05)-biasmean) * dn2msb
print,'1DN 1 sec L1.0 (MSB) mean = ',mean(data_1DN_10)
print,'128DN 1 sec L0.5 (DN) mean = ',mean(data_128DN_05),$
  ' -biasmean * dn2msb = ',(mean(data_128DN_05) -biasmean) * dn2msb
print,'128DN 1 sec L1.0 (MSB) mean = ',mean(data_128DN_10)
print,'1DN 100 sec L1.0 (MSB) mean = ',mean(data100sec),$
  ' * 100 = ',100*mean(data100sec)
print,'So, in summary, these five should all be roughly equal:'
print,'dn2msb, (1DN_05-biasmean)*dn2msb, 1DN_10, 128DN_10/128, data100sec*100'
print,mean(dn2msb*data_flat),$
  (mean(data_1DN_05)-biasmean) * dn2msb,$
  mean(data_1DN_10),$
  mean(data_128DN_10)/128,$
  mean(data100sec)*100

  if (plot eq 0) then stop

; #############################
; # Plot the data and some ratios
  !p.multi = [0, 2, 4, 0, 1]
  plot,data_1DN_05c,title='1DN, level 0.5'
  plot,data_1DN_10c,title='1DN, level 1.0'
  plot,data_1DN_05c/data_1DN_10c,title='1DN, ratio 0.5/1.0'
  plot,data_flatc,title='flat field cal file'

  plot,data_128DN_05c,title='128DN, level 0.5'
  plot,data_128DN_10c,title='128DN, level 1.0'
  plot,data_128DN_05c/data_128DN_10c,title='128DN, ratio 0.5/1.0'
  plot,data_flatc,title='flat field cal file'
  !p.multi=0
; # and save it
  write_gif,'1,128DN:L0.5,1.0.gif',tvrd()

; # plot the images of the data
  datums=[[[data_1DN_05]],[[data_1DN_10]],[[data_128DN_05]],[[data_128DN_10]],[[data_flat]]]
  tv_multi,datums,title='1, 1, 128, 128',/log,res=0.1

  datums=[[[data_1DN_05c]],[[data_1DN_10c]],[[data_128DN_05c]],[[data_128DN_10c]],[[data_flatc]]]
  tv_multi,datums,title='1, 1, 128, 128, congrid128',/log,res=0.1

; # Now plot our 1 and 100 sec versus the flat field

  data100sec=sccreadfits(fname1DN100secL10)
  !p.multi=[0,3]
  plot,alog10(data100sec[*,1023]),title='100 sec L1 B/B0',yrange=[-13.5,-10]
  plot,alog10(data_1DN_10[*,1023]),title='1 sec L1 B/B0',yrange=[-13.5,-10]
  plot,alog10(data_msb[*,1023]),title='MSB*flat',yrange=[-13.5,-10]
  !p.multi=0
  write_gif,'1vs100sec.gif',tvrd()

; -10.5 to -5.5
  plot,alog10(data100sec[*,1023]),yrange=[-13.5,-10],$
    title='L1.0 1DN 1 sec, 100sec, vs MSB'
  oplot,alog10(data_1DN_10[*,1023])
  oplot,alog10(data_msb[*,1023])
  write_gif,'1vs100secalt.gif',tvrd()

  !p.multi=[0,2,2]
  plot,alog10( (data_1DN_05[*,1023]-biasmean)*dn2msb*data_flat[*,1023] ),$
    yrange=[-11.5,-10],$
    title='B/B0: 1DN L0.5 - biasmean * dn2msb * flat'
  plot,alog10(data_1DN_10[*,1023]),yrange=[-11.5,-10],$
    title='B/B0: 1DN L1.0'
  plot,alog10(data_128DN_10[*,1023]/128.),yrange=[-11.5,-10],$
    title='B/B0: 128DN L1.0 /128'
  plot,alog10(data100sec[*,1023]*100.),yrange=[-11.5,-10],$
    title='B/B0: 1DN 100 sec L1.0 * 100'
  !p.multi=0
  write_gif,'B0verify1.gif',tvrd()

  plot,alog10( (data_1DN_05[*,1023]-biasmean)*dn2msb*data_flat[*,1023] ),$
    yrange=[-11.7,-10],title='4 overlays of MSB'
  oplot,alog10(data_1DN_10[*,1023])
  oplot,alog10(data_128DN_10[*,1023]/128.)
  oplot,alog10(data100sec[*,1023]*100.)
  write_gif,'B0verify2.gif',tvrd()