function dnpersec,Is,temps,wave

; CSee  27 Feb 2007, rev 07May07
; This program creates SA data numbers (DN/sec) from Lyn's Titan solar intensity model.
; This version incorporates the red SA channels.
; The model should be interpolated spatially to agree with the pointing angle of the SA pixels.
; This is the function version which pass in:
; 1) Is(col,row,Nwaves,chain) - Solar intensity for each SA pixel verses wavelength. fltarr(6,50,Nwave,4)
; 2) Temps(2) - CCD and Optics Temperatures (float)
; 3) wave - An array specifing the wavelengths used in Lyn's model. (and range of integration)
; And returns DN(col,row,chain) - The DN/sec in the SA for each pixel and chain (see note below)
; Note: chain 0=bv, 1=bh, 2=rv, 3=rh

;_____________________________________________
;Approach....

;Equations...
;1) DN/sec = Abs_Resp * Integral_Lambda(Pol_Factor * Model_Intensity * RSR) d_Lambda
;2)  Abs_Resp = Mean_Abs_Resp_polynomial * Per_Pixel_Polynomial
;3)   Abs_Resp = (A0 + A1*TO + A2*TO^2 + A3*TO^3 + A4*TO^4 + A5*TO^5 + A6*TO^6) * $
;                (MA0 + MA1*T + MA2*T^2 + MA3*T^3) ;per row & column
;      where T is the detector temperature and TO is the optics temperature.
;4)  RSR = RSR_filter * Bump * RSR_NFM  ; NFM is made up of temperature dependent
;                                       ; and non-temp dependent terms
;5)   RSR(lambda,row,column,T) = RSR_filter * Bump * [(M0 + M1*Lambda + M2*Lambda^2 + M3*Lambda^3) + $
;           (224.7-T) * (a + b*Lambda + c*Lambda^2 + d*Lambda^3)] ; for each column
;6)   For the Blue SA, Bump = 1.
;     For the Red SA, Bump = Bump_Lambda * (nfm0 + nfm1*Lambda + nfm2*Lambda^2)
;7)    Bump_Lambda = 1+(0.4/(1+[2*(854-Lambda)/15]^2)
;8)  Polorization_Factor = 0.5  ;assumption for non-polorized light near the sun

;Steps...
; Do one channel at a time (blue_vertical, blue_horizontal, etc)
;
;1) Accepts the Intensities, I(col,row,lambda,chain); Tenperatures, (CCD,Optics);  and model
;    wavelengths, wave(Nwave) as pass parameters
;
;2) Calculate Absolute Responsivity...
; a) Read in Mean_Abs_Resp_polynomial cofficients (A's) from mean_abs_resp_coefs.prn
; b) Calculate Mean_Abs_Resp_Polynomial (marp) at optics temperature TO
; c) Read in Per_Pixel_Polynomial coeficients [MAx(col,row)] from file
;     (i.e. abs_resp_ppp_blue_horiz.txt)
; d) Calculate Absolute_Responsivity [Abs_Resp(col,row)] at the model temperatures, T & TO
;     (equation 3 above)
;
;3) Calculate Relative Spectral Response...
; a) For red SA calculate bump(lambda) from equations 6 & 7.
; b) Read in filter data, FilterTable(row,lambda) from file (i.e. filter_curves_blue_vert.prn)
; c) Interpolate to Intensity and SA row grids [filt(row,lambda)]
; d) Read-in M coeficients, Mx(column), from file (i.e. rsr_nfm_coef_blue_horiz.prn)
; e) Read-in temperature cofficients, B(i) from file (rsr_nfm_temp_coefs.prn) a=B0, b=B1, etc.
; f) Calculate RSR using equations 4 & 5 above
;
;4) Create the integral for each pixel...
; a) for each column, col=0,5 do begin
; b)  for each row, row=0,49 do begin
; c)   sum over wavlength range: for i=lambda_1,lambda_n do begin
; d)    Int(col,row)=Int(col,row) + I(col,row,lambda)/2 * Bump(lambda) * RSR(column,lambda) * $
;        [lambda(i)-lambda(i-1)]
; e) end lambda, end row, end col

;5) Calculate DN/sec for each pixel...
; a) DN(col,row)=Abs_Resp(col,row) * Int(col,row)



start:



;___________________________________________________
;Initialize...
;;print,'reading-in parameters...'
;;wait,0.5

DN=dblarr(6,50,4)  ;DN_per_second(col,row,chain)
pol_fact = 0.5     ; assumption that 1/2 the signal is in each chain for unpolarized light.
line=''            ; variable for reading from files.
T=temps(0)	    ;CCD Temp
TO=temps(1)	    ;Optics Temp
s=size(IS)
Nwave=s(3)	    ;number of wavelengths
ch1=0           ;first channel (i.e. Ch1=0 is blue vertical)
ch2=3           ;last channel to process (Ch2=3 is red horizontal)

if s(4) eq 2 then begin   ;i.e. only 2 channels (blue or red only)
 if wave(0) lt 650 then ch1=0 else ch1=2   ;blue vs red
 ch2=ch1+1
endif


for chain=ch1,ch2 do begin	;0=bv, 1=bh, 2=rv, 3=rh

;Files...
dir=                       'D:\Documents\disr\SA_simulator\'    ;Working Directory
mean_abs_resp_coefs_file=  'mean_abs_resp_coefs.prn'            ;3) equation
rsr_nfm_temp_coefs_file=   'rsr_nfm_temp_coefs.prn'             ;5)

if chain eq 0 then begin
 per_pixel_poly_coefs_file= 'abs_resp_ppp_blue_vert.txt'        ;3)
 filter_curve_file=         'filter_curves_blue_vert.prn'       ;5)
 rsr_nfm_file=              'rsr_nfm_coef_blue_vert.prn'        ;5)
endif

if chain eq 1 then begin
 per_pixel_poly_coefs_file= 'abs_resp_ppp_blue_horiz.txt'      ;3)
 filter_curve_file=         'filter_curves_blue_horiz.prn'     ;5)
 rsr_nfm_file=              'rsr_nfm_coef_blue_horiz.prn'      ;5)
endif

if chain eq 2 then begin
 per_pixel_poly_coefs_file= 'abs_resp_ppp_red_vert.txt'         ;3)
 filter_curve_file=         'filter_curves_red_vert.prn'        ;5)
 rsr_nfm_file=              'rsr_nfm_coef_red_vert.prn'         ;5)
endif

if chain eq 3 then begin
 per_pixel_poly_coefs_file= 'abs_resp_ppp_red_horiz.txt'       ;3)
 filter_curve_file=         'filter_curves_red_horiz.prn'      ;5)
 rsr_nfm_file=              'rsr_nfm_coef_red_horiz.prn'       ;5)
endif

bump_file=                 'red_bump_coefs.prn'	                ;6) only for red SA






;___________________________________________________
;Calculate Absolute Responsivity...

;read-in mean_abs_resp_coefs...
close,3
openr,3,dir+mean_abs_resp_coefs_file
Aa=dblarr(7)	;array containing the AbsResp coefs

case chain of
 1: index=1		;bh
 0: index=2		;bv
 3: index=3		;rh
 2: index=4		;rv
endcase

readf,3,line    ;header
for i=0,6 do begin
 readf,3,line
 Aa(i)=double(strmid(line,20*index,20))
endfor
close,3

;Calculate the Mean Absolute Responsivity at this optics temperature (TO)...
marp=(Aa(0) + Aa(1)*TO + Aa(2)*TO^2 + Aa(3)*TO^3 + Aa(4)*TO^4+ Aa(5)*TO^5 + Aa(6)*TO^6)

;read-in the per-pixel abs resp coefs...
close,4
openr,4,dir+per_pixel_poly_coefs_file
Ma=dblarr(4,6,50)   ; per pixel coeficients M0,M1,M2,M3 {Ma(coef,col,row)}
For i=0,3 do begin  ; for each coefficient...
 readf,4,line       ; 1st header
 readf,4,line       ; 2nd header
 for row=0,49 do begin
  readf,4,line
  for col=0,5 do Ma(i,col,row)=double(strmid(line,4+21*col,20))
 endfor
endfor
close,4

;calculate the absolute responsivity...
Abs_Resp=dblarr(6,50)  ; for each pixel, Abs_Resp(col,row)
for col=0,5 do begin
 for row=0,49 do $
  Abs_Resp(col,row)= marp * (Ma(0,col,row) + Ma(1,col,row)*T + Ma(2,col,row)*T^2 + Ma(3,col,row)*T^3)
endfor


;___________________________________________________
;Calculate Relative Spectral Response...

;calculate bump...
bump=dblarr(nwave)
if chain eq 0 or chain eq 1 then bump(*)=1.0 $	;always 1 for blue SA
else begin	                                ;for the red SA...
 Bump_Lambda = 1+(0.4/(1+(2*(854-wave)/15)^2))
 ;read-in bump_coefs...
  close,8
  openr,8,dir+bump_file
  readf,8,line    ;header
  readf,8,line    ;nfm0
  if chain eq 3 then nfm0=double(strmid(line,10,20)) else nfm0=double(strmid(line,30,20))
  readf,8,line    ;nfm1
  if chain eq 3 then nfm1=double(strmid(line,10,20)) else nfm1=double(strmid(line,30,20))
  readf,8,line    ;nfm2
  if chain eq 3 then nfm2=double(strmid(line,10,20)) else nfm2=double(strmid(line,30,20))
  close,8
 Bump = Bump_Lambda * (nfm0 + nfm1*wave + nfm2*wave^2)
endelse


;read in filter curves...
if chain eq 0 or chain eq 1 then begin   ;blue SA
 nCurves=13                               ; 13 samples in SA row dimension
 nFWaves=37                               ; the number of wavelengths in the filter table
 rowf=fltarr(nCurves)                     ; an array to hold the filter's equiv. row number
 fwave=intarr(nfwaves)                    ; an array to hold the tables wavelengths
 FilterTable=fltarr(nCurves,nFWaves)      ; filter table data
 temp=fltarr(nCurves,nwave)               ; temporary file for interpolation to solar I grid
 Filt=fltarr(50,nwave)       ; the filter transmission for each row vs. wavelength
endif else begin                         ;red SA
 nCurves=13                               ; 13 samples in SA row dimension
 nFWaves=43                               ; the number of wavelengths in the filter table
 rowf=fltarr(nCurves)                     ; an array to hold the filter's equiv. row number
 fwave=intarr(nfwaves)                    ; an array to hold the tables wavelengths
 FilterTable=fltarr(nCurves,nFWaves)      ; filter table data
 temp=fltarr(nCurves,nwave)               ; temporary file for interpolation to solar I grid
 Filt=fltarr(50,nwave)       ; the filter transmission for each row vs. wavelength
endelse

close,5
openr,5,dir+filter_curve_file
readf,5,line   ;list of row numbers
for i=0,nCurves-1 do rowf(i)= float(strmid(line,10+10*i,10))
 check=[0.5,5.0,10.0,15.0,20.0,25.0,30.0,35.0,40.0,43.0,45.0,48.5,50.0]
 if where(rowf ne check) ne -1 then print,'*** THE FILTER CURVE FILES APPEARS TO BE WRONG ***'
readf,5,line   ; a header
for i=0,nFWaves-1 do begin
 readf,5,line
 fwave(i)=float(strmid(line,0,10))
 for j=0,ncurves-1 do FilterTable(j,i)=float(strmid(line,10+10*j))
endfor
close,5

;interpolate filter data in wavelength to match Solar Intensity wavelenght grid.
for k=0,ncurves-1 do begin
 temp(k,*)=interpol(filterTable(k,*),fwave,wave)
endfor

;create filter data for each row...
row=findgen(50)
for m=0,nwave-1 do filt(*,m)=interpol(temp(*,m),rowf,row)

;read-in M coefs...
Mx=dblarr(4,6)   ; RSR M coefs, M0,M1,M2,M3 {Mx(i,column)}
close,6
openr,6,dir+rsr_nfm_file
readf,6,line		;header
for col=0,5 do begin
 readf,6,line
 for i=0,3 do Mx(i,col)=double(strmid(line,20+i*20,20))
endfor
close,6

;read-in RSR NFM Temperature, B coefs (a=B0, b=B1, etc)
Bx=dblarr(4)	; a,b,c & d
close,7
openr,7,dir+rsr_nfm_temp_coefs_file
readf,7,line   ;header

case chain of
 1: index=1    ;bh
 0: index=2    ;bv
 3: index=3    ;rh
 2: index=4    ;rv
endcase

for i=1,index do readf,7,line
for j=0,3 do Bx(j)=double(strmid(line,20+20*j,20))

if chain eq 3 then begin
 B0=dblarr(4)	;special coef for red vertical column 0
 B1=dblarr(4)	;special coef for red vertical column 1
 readf,7,line   ;red vertical
 readf,7,line   ;red horizontal_col_0
  for j=0,3 do B0(j)=double(strmid(line,20+20*j,20))
 readf,7,line   ;red horizontal_col_1
  for j=0,3 do B1(j)=double(strmid(line,20+20*j,20))
endif
close,7

;calculate RSR Non-Filter Model (RSR_NFM) as in equation 5...
rsr_nfm=dblarr(6,nwave)   ;RSR_NFM(col,lambda)
for col=0,5 do rsr_nfm(col,*)=(Mx(0,col) + Mx(1,col)*wave + Mx(2,col)*wave^2 + Mx(3,col)*wave^3) + $
 (224.7-T) * [Bx(0) + Bx(1)*wave + Bx(2)*wave^2 + Bx(3)*wave^3]

;calculate RSR at this temperature (eqn 4)...
rsr=dblarr(6,50,nwave)   ;Relative_Spectral_Response(Column,Row,Lambda)
for col=0,5 do begin
 for row=0,49 do rsr(col,row,*)=filt(row,*) * bump(*) * rsr_nfm(col,*)
endfor


;___________________________________________________
;Create Integral of Solar_Intensity * Pol_Factor * RSR dLambda...
;;print,'Integrating...'
;;wait,0.5
Int=dblarr(6,50)	;Integral(col,row)
for col=0,5 do begin
 for row=0,49 do begin
  for k=1,nwave-1 do Int(col,row)=Int(col,row) + $
   pol_fact * 0.5*[Is(col,row,k-1,chain)*rsr(col,row,k-1)+ $
   Is(col,row,k,chain)*rsr(col,row,k)] * bump(k) * [wave(k)-wave(k-1)]/1000  ;polygon integral
 endfor  ;row
endfor  ;col


;___________________________________________________
;Calculate DN/sec (eqn 1)...
;;print,'Calculating DN/sec...'
;;wait,0.5


DN(*,*,chain)=abs_resp * Int

endfor   ;chain

exit:
;Print,string(10b),'Done!  :)'
;stop
close,/all
return,DN
end