;+
; Project     : SOHO - CDS     
;                   
; Name        : VDS_DUMMY
;               
; Purpose     : Produce dummy VDS data in the form of a (1024,256) array.
;               
; Explanation : Takes past solar EUV data, folds in CDS performance
;               and solar parameters to mimic 1024x256 portion of
;               VDS image. Includes -
;                    (i)    solar lines in correct location on image
;                    (ii)   intensities calculated from past data
;                           folding in CDS efficiencies
;                    (iii)  effects of telescope point spread function (3" FWHM)
;                    (iv)   variable slit use
;                    (v)    several window options
;                    (vi)   dummy solar structure along slit
;                    (vii)  background intensities
;                    (viii) statistical noise
;                    (ix)   Doppler shifted ``events''
;
;               Neglects pixel sampling effects (assumed small) and
;               grating scatter (assumed small w.r.t. slit image
;               spread).
;
;               If requested, the output array is written to a disk file
;               with an unformatted write.  In order to retrieve the data
;               use the following commands:
;               IDL> a = intarr(1024,256)
;               IDL> openr,lun,filename,/get_lun
;               IDL> readu,lun,a
;               IDL> free_lun,lun
;          
; Use         : IDL> vds_dummy [,slit=slit,expose=expose etc.]
;               IDL> vds_dummy, slit=4,expose=10,wlist=1,file='vds1.dat',
;                               image=vds1
;               would return the image in variable vds1 and store in a file
;               called vds1.dat
;
; Inputs      : Slit width, Accumulation time and Window option
;               are prompted for if not supplied as keyword parameters.
;               
; Opt. Inputs : None
;               
; Outputs     : None
;               
; Opt. Outputs: See Image keyword
;               
; Keywords    : Slit         The slit size (2,4,8 or 90 arcsec).
;               Expose       The exposure time in seconds.
;               Wlist        The line list reference number.
;               Win_width    Individual line window width in pixels.
;               Image        Name of data array returned.
;               File         Name of data file created (default image.dat).
;               No_file      Inhibit saving of data on disk.
;               Zone         Defines type of source model.  The only line
;                            intensity data available at present are for
;                            QUIET regions (zone=4) and ACTIVE regions (zone=5)
;
; Calls       :
;               
; Restrictions: 1025 blocks required for output file if any.
;               Requires solar data in files ni1.dat, ni2.dat and
;               pixel address information in pixels.dat
;               
; Side effects: None
;               
; Category    : Planning, VDS, Model
;               
; Prev. Hist. : None
;
; Written     : R A Harrison, RAL, March 1993
;               
; Modified    : Version 2, William Thompson, GSFC, 17 May 1993.
;                     Corrected algorithm used for spectral dispersion.
;                     Added dispersion within VDS detector.
;		Version 3, William Thompson, GSFC, 22 September 1993.
;			Fixed typo in dispersion calculation.
;
;               Version 4, CDP, Update format of data files to handle 6 zones.
;                                       27-March-95
;               Version 5, CDP, Manual input of zone if not on command line
;                                       22-Aug-95
;               Version 6, trap invalid ZONE IDs. CDP, 27-Nov-95
;
;               Version 7, Now includes the O IV 629 A line P.B., 18-Feb-96
;
;               Version 8, Fix bug in manual zone ID input. CDP, 20-Feb-96
;
;               Version 9, Bump up (down?) spectral resolution to match
;                          reality.  CDP, 27-Apr-96
;  
; Version     : Version 9, 7 April 1996 
;-            

pro vds_dummy,slit=slit,expose=expose,wlist=wlist,file=file,$
              image=image,win_width=win_width,no_file=no_file,$
              debug=debug,zone=zone

;
;  Trap i/o errors
;
on_ioerror,finish
;
;  File contains pixel number, upper spectrum, lower spectrum wavelengths
;
band1=fltarr(1024,2)       ;1024 pixels 310-380A band, 0=wavelength, 1=intensity
band2=fltarr(1024,2)       ;ditto 517-633A band

band1(*,0)=pix2wave('N1',indgen(1024))
band2(*,0)=pix2wave('N2',indgen(1024))

;
; Request slit size if not set on command line
;
good_slit = [2,4,8,90]
if not keyword_set(slit) then begin
   read,' Please input slit width choice - 2, 4, 8, 90 arcsec .....>',slit
endif
;
;  Check valid response
;
ok = where(slit eq good_slit)
while ok(0) lt 0 do begin
   read,' Please input slit width choice - 2, 4, 8, 90 arcsec .....>',slit
   ok = where(slit eq good_slit)
endwhile

slit=slit/2.0              ;slit width in units of pixels
if keyword_set(debug) then stop,0

;
;  is zone specified?
;
if not keyword_set(zone) then begin
   read,'Enter Zone ID (4=Quiet, 5=Active): ', use_zone
   if use_zone ne 4 and use_zone ne 5 then begin
      print,'Dumbo!'
      return
   endif
endif else begin
   use_zone = zone
endelse

;
;  for now the maximum allowed zone ID is 5
;
if use_zone gt 5 then begin
   print,'Invalid zone ID for vds_dummy.  Set to <QUIET>.'
   use_zone = 4
endif 

;
;  Fix telescope psf
;
tel=sqrt(2.)*0.637       ;sigma of telescope function (i.e. FWHM=3.0 pixels)
                     ;convolved with detector resolution (also FWHM=3.0 pixels)

grt=sqrt(0.667^2 + 0.637^2)
                     ;sigma of line width (assumed to be 110 mA from SERTS)
                     ;convolved with detector resolution

;
;  increase FWHM to match reality
;
grt = grt * 1.9

;
; Calculate Gaussian profiles for spectral(x) and spatial(y) dimensions
; store spectral profile in array 101 which is barely enough for slit 90 
; but overkill for others and spatial profile in array 11
;
xprofile=fltarr(101)       ;Gaussian profile from slit
yprofile=fltarr(11)       ;Guassian profile from telescope

;
;  In the X direction convolve the Gaussian line profile (which already
;  contains the smearing due to the detector) with the shape of the slit.
;  Note: For even slit widths (in pixels), this will shift the centroid of the
;  line 1/2 pixel to the right.
;
xi = indgen(101)-50
xprofile = convol( exp(-xi^2/(2.*grt*grt)), replicate(1.,slit))
;
;  In the Y direction simply use the combined Gaussian profile due to the
;  telescope and detector.
;
xi = indgen(11)-5
yprofile = exp(-xi^2/(2.*tel*tel))

;
;  Normalize intensity
;
xprofile = xprofile / (sqrt(2.*!pi) * grt)
yprofile = yprofile / (sqrt(2.*!pi) * tel)
if keyword_set(debug) then stop,1

;
;  Read the line intensities in upper spectrum
;
openr,lun_ni1,concat_dir('$CDS_PLAN_TECH','ni1.dat'),/get_lun   ;310-380A intensities
count = 0
wave1 = fltarr(500) & inten1 = fltarr(500) & datum=fltarr(7)

while not eof(lun_ni1) do begin
   readf,lun_ni1,datum
   wave1(count)  = datum(0)
   inten1(count) = datum(use_zone)
   count = count + 1
endwhile
free_lun,lun_ni1
;
;  Count number found and trim arrays
;
count = where(wave1 gt 0.0)
nline1 = n_elements(count)
wave1  =  wave1(0:nline1-1)
inten1 = inten1(0:nline1-1)
if keyword_set(debug) then stop,2
;
;  ... and the lower spectrum
;
openr,lun_ni2,concat_dir('$CDS_PLAN_TECH','ni2.dat'),/get_lun ;517-633A intensities
count = 0
wave2 = fltarr(500) & inten2 = fltarr(500)
while not eof(lun_ni2) do begin
   readf,lun_ni1,datum
   wave2(count)  = datum(0)
   inten2(count) = datum(use_zone)
   count = count + 1
endwhile
free_lun,lun_ni2

count = where(wave2 gt 0.0)
nline2 = n_elements(count)
wave2  =  wave2(0:nline2-1)
inten2 = inten2(0:nline2-1)
if keyword_set(debug) then stop,3
;
;  Find where line falls in wavelength array and convolve the xprofile
;  with the intensity
;
;  First the upper spectrum
;
for i=0,nline1-1 do begin
   n = max(where(band1(*,0) le wave1(i))) 
   if n ge 50 and n le 973 then begin
      band1(n-50:n+50,1) = band1(n-50:n+50,1) + xprofile * inten1(i)
   endif
endfor

;
;  ... and then the lower
;      Changed slightly here to include the O IV line at 629 A (pixel 994)
;
for i=0,nline2-1 do begin
   n = max(where(band2(*,0) le wave2(i)))
   if n ge 50 and n le 996 then begin     
      band2(n-50:n+20,1) = band2(n-50:n+20,1) + xprofile * inten2(i)
   endif
endfor
if keyword_set(debug) then stop,3.5
;
;  Get the accumulation time if it was not set on the command line
;
if not keyword_set(expose) then begin
   read,' Input Accumulation Time (Sec) .........>',expose
endif
;
;  Bump up counts accordingly
;
;plot,band2,xr=[1000,2050],xstyle=1,yr=[0,50]
band1(*,1)=band1(*,1)*expose
band2(*,1)=band2(*,1)*expose


;
;  Get line list required if not specified on command line
;
if not keyword_set(wlist) then begin
   print,' Input Window/Line List - Options are:'
   print,'     1   all data'
   print,'     2   list L12 (Blue Book) with 21 pixels across each line'
   print,'     3   two windows - 345-360A and 550-565A'
   read,'      .........>',wlist
endif
good_wlist = [1,2,3]
;
;  Check valid response
;
ok = where(wlist eq good_wlist)
while ok(0) lt 0 do begin
;  bell
   read,' Please input valid option: ',wlist
   ok = where(wlist eq good_wlist)
endwhile
;
;  Set up mask according to line list option chosen
;
;  Assume window width of 21 pixels unless specified as keyword
;

mask1=bytarr(1024)       ;upper window mask
mask2=bytarr(1024)       ;lower window mask

if not keyword_set(win_width) then win_width = 21
case wlist of 
   1: begin 
         mask1(*) = 1 
         mask2(*) = 1
      end
   2: begin
         b = [355,447,488,561,619,821]         ; upper spectrum
         b1 = b
         for i=1,win_width do b = [b,b1+i]
         mask1(b) = 1
         b = [21,27,310,535,564,732,754,762,910]  ; lower spectrum
         b1 = b
         for i=1,win_width do b = [b,b1+i]
         mask2(b) = 1
      end
   3: begin
         mask1(502:715) = 1
         mask2(282:408) = 1
      end
   else: print,'Invalid line list option:'
endcase
if keyword_set(debug) then stop,4
;
;  Add an assumed background of 1 c/s and mask off the data
;
band1(*,1) = (band1(*,1)+expose) * mask1
band2(*,1) = (band2(*,1)+expose) * mask2

;
; Input Dummy Solar Structure Along Slit
;
structure = fltarr(120) + 1.0
for i=0,119,10 do begin
   structure(i)=structure(i)*1.5
endfor
for i=25,75,3 do begin
   structure(i)=structure(i)*2.2
endfor
for i=0,119,17 do begin
   structure(i)=structure(i)*1.7
endfor
for i=0,119,3 do begin
   structure(i)=structure(i)*1.1
endfor
if keyword_set(debug) then stop,5
;
;  Put in Gaussian for Telescope Point Spread Function
;
structure = convol(structure,yprofile)
;
;  Patch up the ends, missed because of width of kernel array yprofile
;
structure(0:4) = structure(5)
structure(115:119) = structure(114)

if keyword_set(debug) then stop,6
;
;  Read into 2-D array
;
image = fltarr(1024,256)
image(0:1023,9:128)   =  band2(*,1) # float(structure) 
image(0:1023,128:247) =  band1(*,1) # float(structure) 

if keyword_set(debug) then stop,7

;
;  Add noise - random distribution, function of intensity
;
a  = (randomu(seed,1024,256) - 0.5) * 2.0
image = image + (a * sqrt(image))
;
;  Force to be positive
;
image = image > 0

if keyword_set(debug) then stop,8

;
;  Add Doppler shifts - 30 randomly selected locations
;
;  5 single pixel positive
;
a = fix(randomu(seed,4) * 255)
image(*,a) = 0.4*image(*,a) + 0.4*shift(image(*,a),1) + $
             0.2*shift(image(*,a),2)
;
;  5 single pixel negative
;
a = fix(randomu(seed,4) * 255)
image(*,a) = 0.5*image(*,a) + 0.3*shift(image(*,a),-1) + $
             0.2*shift(image(*,a),-2)
;
; 10, 3 pixel positive
;
a = fix(randomu(seed,10) * 253)
image(*,a)   = 0.7*image(*,a)   + 0.2*shift(image(*,a),1)   + $
                                  0.1*shift(image(*,a),2)
image(*,a+1) = 0.6*image(*,a+1) + 0.3*shift(image(*,a+1),1) + $
                                  0.2*shift(image(*,a+1),2)
image(*,a+2) = 0.6*image(*,a+2) + 0.2*shift(image(*,a+2),1) + $
                                  0.1*shift(image(*,a+2),2)
;
; 10, 3 pixel negative
;
a = fix(randomu(seed,10) * 253)
image(*,a)   = 0.5*image(*,a)   + 0.2*shift(image(*,a),-1)   + $
                                  0.2*shift(image(*,a),-2)
image(*,a+1) = 0.6*image(*,a+1) + 0.3*shift(image(*,a+1),-1) + $
                                  0.2*shift(image(*,a+1),-2)
image(*,a+2) = 0.5*image(*,a+2) + 0.3*shift(image(*,a+2),-1) + $
                                  0.2*shift(image(*,a+2),-2)

;
;  Fix and trim image array if needed elsewhere
;
if keyword_set(image) or (not keyword_set(no_file)) then begin
   image = fix(image < 4096)
endif
;
;  Create output file if needed.
;
if keyword_set(debug) then stop,9
if not keyword_set(no_file) then begin
   if keyword_set(file) then fname = file else fname = 'image.dat'
   print,'Output file name is '+fname
   openw,lun_out,fname,/get_lun
   writeu,lun_out,image
   free_lun,lun_out
endif

print,' VDS_Dummy - Normal End'
return

;
;  error return
;
finish:
;bell
print,' VDS_Dummy - error return'
print,'Probably trouble reading the input data files:'
print,'ni1.dat and ni2.dat'
print,'which should be in directory $CDS_PLAN_TECH'
print,' '

end