pro OptimizeV9,alt,ModNum

alt=80
ALT=88
alt=100
alt=60.4481	;test case
color='blue'    ;'red', 'blue', or 'both' <- Set SA Model color here

q=systime(1)   ;how long does it take?

;CSee 26Feb07, revised 03July2007

;This program finds the optimum match between Lyn's model data and the Titan
; descent data.  It processes all models in the Model_directory (or a single
; model if designated by ModNum) and minimizes the standard deviation of the
; difference between calculated data numbers (using SA_sim.pro) and the actual
; Titan descent data (dark subtracted), in a 3 dimensional parameter space of
; Azimuth (Az), Tip toward the sun (T1), and Tip perpendicular to T1 (T2).
;
;The sign convention is Az positive CCW (from top) relative to the sun, T1
; positive for the + probe spin axis tipping away from the sun, and T2 positive
; is the probe spin axis tipping southwestward (RHR).  Note: a positive T1 or
; and/or T2 results in a westward tip, i.e. the opposite sign from Erich's
; convention.
;
;Input parameters...
; Alt=target altitude for analysis in KM (floating)
; ModNum=Sequential number for running just an individual number (integer)
;
;Approach...
;1) Find the lamps-off dataset closest to the target altitude and with azimuth
;    close to the sun
;2) Define initial point using Az from header and Ti=T2=0
;3) Create a grid of values using a range and steps defined below.
;4) Run SA simulator at each point and calculate RMS deviation (data to simulator)
;    at each pixel.  Nested loops, T2 inner, T1 second, AZ outer.
;5) Keep the points bounding the minimum (RMS deviation) as endpoints and repeat 3
;    with a tighter grid.
;6) Save Closest Result, Az, T1, T2, Min_RMS_dev, DN/sec, Difference Array
;7) Repeat for next model.


;____________________
;Initilaization...

; Y=//laplace/data/ldoose
; Z=//Cassini/descent
;model_dir='Y:\Titan_ULVS_models\Model_777'  ;directory containing blue model
;model_dir='Y:\Titan_ULVS_models\Model_853'  ;directory containing red model
;model_dir='Y:\Titan_ULVS_models\Model_736\Red'  ;directory containing red model (numg=11)
;model_dir='Y:\Titan_ULVS_models\Model_860'  ;directory containing red model
model_dir='Y:\Titan_ULVS_models\Model_861'  ;directory containing blue model
data_dir='Z:\Log\stream_524a_alt_az_5-8-2006\DB\Solar'
Log='Z:\Log\stream_524a_alt_az_5-8-2006'
dir='D:\Documents\disr\SA_simulator'	;working directory (files + output)


;Which SA channel, blue or red...
  ch1=0 & ch2=3	;ch1 is the first chain (blue vert) and ch2 is the last (red horiz)
  if color ne 'both' then if color eq 'blue' then ch2=1 else ch1=2

elapsed=fltarr(6)  ;keep track of elapsed time in each section
oldt=systime(1)

AzRangeToSun=40. ; Restriction on azimuth angle to the sun when selecting data.




common saved,model_name_old,alt_old,tables_read,wavl,solar_zen_angle, $
  theta,phi,idnc,zeniths_blue_vertical,azimuths_blue_vertical_init, $
  zeniths_blue_horizontal,azimuths_blue_horizontal_init

; These common block variables must be initialized
model_name_old=''
alt_old=-1.0
tables_read=0
solar_az=336.0
solar_zen_angle=36.0
sza=solar_zen_angle
pi=3.141592654d
esa=24.*pi/180	;East Solar Angle (Rads), angle of the sun relative to east (114-90)




;____________________
print,'Finding Descent Data...'
wait,0.5

 ;find which Descent data to use (closest to target altitude)...
 files=file_search(data_dir,'*MMX*',count=nfiles)
 array_sort,files	;sorts 'files' into assending order
 seqno=0       ; sequence number of current data file
 altdata=0.0   ; The altitude of current data file
 for fileno=0,nfiles-1 do begin
  d_read,files(fileno),h,p
  ;skip data with lamps on or too far from sun...
   az=d_value(h,101)    ;h101=actual_azimuth
   if az le 174 then AngleToSun=az+6 else AngleToSun=abs(6-(360-AZ))
   if d_value(h,102) ne '0000' or AngleToSun gt AzRangeToSun then goto,nextfile
   if d_value(h,6) ne 24 then goto,nextfile   ;columns = 24
   case color of
    'blue': if max(p(0:11,*)) gt 4000 then goto,nextfile       ;saturated data
    'red' : if max(p(12:23)) gt 4000 then goto,nextfile       ;saturated data
    'both': if max(p) gt 4000 then goto,nextfile       ;saturated data
   endcase
  curralt=d_value(h,98)/1000.
  currseqno=d_value(h,87)

  ;start from top of descent and work down
  if curralt gt alt then begin
   high=curralt     ; above target altitude
   hsqn=currseqno
   hfile=fileno
   goto,nextfile
  endif else begin
    low=curralt     ; below target altitude
    lsqn=currseqno
    lfile=fileno
    goto,findclosest  ;quit after first dataset below target altitude
  endelse
  nextfile:
 endfor	  ;loop on files

 findclosest:
 if high-alt gt alt-low then begin   ;i.e. the lower one is closer to target
  alt=low      ;the sellected altitude
  seqno=lsqn   ;sequence number of the sellected dataset
  fileno=lfile ;index of the lower file
  endif else begin
   alt=high
   seqno=hsqn
   fileno=hfile   ;index of the higher file
  endelse

;Now that we have the right dataset lets get some details...
d_read,files(fileno),h,p
 cols=d_value(h,6)         ;h6=cols ;24
 rows=d_value(h,7)         ;h7=rows ;50
 Tccd=d_value(h,22)        ;h22=ccd_chip
 Etime=d_value(h,36)       ;h36=exptime
 Mtime=d_value(h,88)       ;h88=mission_time_sec
 cycleno=d_value(h,89)     ;h89=cycle_no
 az=d_value(h,101)         ;h101=actual_azimuth
 Toptics=otemp(Mtime,Log)  ;returns the optics temperature from HK data

 if alt ne d_value(h,98)/1000. then stop          ;h98=altitude_meteres
 if seqno ne d_value(h,87) then stop        ;h87=seq#
 if d_value(h,102) ne'0000' then stop  ;h102=lamp_state

temps=[Tccd,Toptics]

;remove dark current
;p=p-8

     ; Remove pixels that for which there is substantial error between the IS data and the SA simulator.
     PixelMaskBV=intarr(6,50)
     pixelmaskbv(*)=1
     ;pixels with errors over 5%...
     pixelmaskbv(0,22)=0
     pixelmaskbv(1,46)=0
     pixelmaskbv(1,48)=0
     pixelmaskbv(4,0)=0
     pixelmaskbv(5,0)=0
     pixelmaskbv(5,46)=0
     pixelmaskbv(5,47)=0
     pixelmaskbv(5,49)=0

     PixelMaskbh=intarr(6,50)
     pixelmaskbh(*)=1
     ;pixels with errors over 5%...
     pixelmaskbh(0,26)=0
     pixelmaskbh(3,47)=0
     pixelmaskbh(3,49)=0
     pixelmaskbh(4,46)=0
     pixelmaskbh(4,47)=0
     pixelmaskbh(4,49)=0
     pixelmaskbh(5,46)=0
     pixelmaskbh(5,47)=0
     pixelmaskbh(5,48)=0
     pixelmaskbh(5,49)=0

     PixelMaskrv=intarr(6,50)
     pixelmaskrv(*)=1

     PixelMaskrh=intarr(6,50)
     pixelmaskrh(*)=1
     pixelmaskrh(0,0:3)=0
     pixelmaskrh(0,16:49)=0
     pixelmaskrh(5,25:49)=0
     pixelmaskrh(1,2)=0
     pixelmaskrh(1,3)=0
     pixelmaskrh(0,10)=0
     pixelmaskrh(0,11)=0
     pixelmaskrh(0,12)=0


    ;Note: there is a polerization inversion in the Blue SA chain...
     ; p(0,5)=Blue Horizontal
     ; p(6,11)=Blue Vertical
     ; p(12,17)=Red Vertical
     ; p(18,23)=Red Horizontal

     pixelmask=[pixelmaskbv,pixelmaskbh,pixelmaskrv,pixelmaskrh]
     pixelmask=pixelmask(6*ch1:6*ch2+5,*)   ;use only the channels in the model (i.e. blue, red or both)

;____________________
;Find model values at same altitude, azimuth, etc...

print,'Looking for Models...'
wait,0.5


;Which models to do...
model=file_search(model_dir,'INTENSI*.OUT',count=Nmodels) ;search for model names
if model eq '' then print,'Models not found'

m1=0
m2=Nmodels-1

;to do just one model...
if n_params() eq 2 then begin
 m1=modnum & m2=modnum
 nmodels=1
endif

for modno=m1,m2 do begin	;loop on models


 ;____________________
 ;********************
 ;Seed state....

 azi=az	;azimuth from descent data
 T1=0.0
 T2=0.0
 model_name=model(modno)
 print,'Model = ',model_name

 az1= 10	 ; negative excursion limit in Azimuth
 az2=+18.         ; positive excursion limit in Azimuth
 azstep=0.5       ; Initial azimuth step size...final is initial devided by 'factor'
 T1a=-20.         ; Limits on T1, T2, Step_size, etc as with Azimuth above...
 T1b=+20.
 T1step=4.
 T2a=-20.
 T2b=+28.
 T2step=4.
 factor1=5.	     ; factor by which azimuth grid is reduce for second pass
 factor2=10.	     ; factor by which t1,t2 grids are reduced for second pass
 ;Hint: azstep should not be greater than 1 and factor should not be greater than 10.
 print,'az,az1,az2,azstep=',az,az1,az2,azstep
 print,'t1a,t1b,t1step=',t1a,t1b,t1step
 print,'t2a,t2b,t2step=',t2a,t2b,t2step
 print,'AZfact,Tfact=',factor1,factor2

 ;capture the low state...
 lowdiffarr=dblarr(12,50)
 lowdn=dblarr(6,50,2)
 lowdev=1e6


 ;____________________
 ;Optimization loop....
 elapsed(0)=elapsed(0)+(systime(1)-oldt)/60.
 oldt=systime(1)

 print,'Optimize...'
 print,systime(0)
 wait,0.5

 pass=1

 Pass:	;begin pass here...

  nazs=fix((az2-az1)/azstep)+1          ;number of azimuth steps
  nt1s=fix((T1b-T1a)/T1step)+1          ;number of T1 steps
  nt2s=fix((T2b-T2a)/T2step)+1          ;number of T2 steps
  dev=fltarr(nmodels,nazs,nt1s,nt2s)	;dev(model,az,T1,T2) - holds rms deviation for each pass
  azarr=fltarr(nazs)
  t1arr=fltarr(nt1s)
  t2arr=fltarr(nt2s)
  minarr=fltarr(nazs)	;to hold the minimum at each azimuth


  azi=az+az1	               ;start with lowest azimuth in range (offset from Erichs estimate)
   if azi ge 360 then azi=azi-360

  for naz=0,nazs-1 do begin    ;for each azimuth
   azarr(naz)=azi

   epsilon=.01
   nt1=0	;index for T1
   for T1i=T1a,T1b+epsilon,T1step do begin
    t1arr(nt1)=t1i

    nt2=0	;index for T2
    for T2i=T2a,T2b+epsilon,T2step do begin
     t2arr(nt2)=t2i

     tip_away_from_sun=T1i
     orthogonal_tip=T2i



     ;;; The components are the tip angle are computed from the nominal probe tip here
     ;;orthogonal_tip=asin(sin(!dtor*probe_tip)*sin(!dtor*solar_az))/!dtor
     ;;tip_away_from_sun=-acos(cos(!dtor*probe_tip)/cos(!dtor*orthogonal_tip))/!dtor

     if azi gt 180 then azimuth=azi-360 else azimuth=azi  ;send routine azimuths in the range +/-180 from sun
     ; First call takes longer because it must read the model
     ;;print,'calculating SA_incident_intensities...'
     ;;wait,0.5

      elapsed(1)=elapsed(1)+(systime(1)-oldt)/60.
      oldt=systime(1)

     SA_incident_intensities,model_name,alt,azimuth,tip_away_from_sun, $
       orthogonal_tip,intensity,dir
     ;if color eq 'blue' then SA_incident_intensities,model_name,alt,azimuth,tip_away_from_sun, $
       ;orthogonal_tip,intensity,dir else SA_incident_intensities_red,model_name,alt,azimuth,tip_away_from_sun, $
       ;orthogonal_tip,intensity,dir   ; for red files with gnum=11

     ;;print,'calculating DN/sec...'
     ;;wait,0.5
      elapsed(2)=elapsed(2)+(systime(1)-oldt)/60.
      oldt=systime(1)

     DN=dnpersecV8(intensity,temps,wavl)	;DN(col,row,chain)

      elapsed(3)=elapsed(3)+(systime(1)-oldt)/60.
      oldt=systime(1)



     ;;print,'calculating differance array...'
     ;;wait,0.5
     diffarr=[DN(*,*,0)-((p(6:11,*)/(etime/1000.))),DN(*,*,1)-((p(0:5,*)/(etime/1000.))), $
              DN(*,*,2)-((p(12:17,*)/(etime/1000.))),DN(*,*,3)-((p(18:23,*)/(etime/1000.)))]

     diffarr=diffarr(6*ch1:6*ch2+5,*)   ;use only the channels in the model (i.e. blue, red or both)
     ;pixelmask=pixelmask(6*ch1:6*ch2+5,*)   ;use only the channels in the model (i.e. blue, red or both)

     diffarr(where(pixelmask eq 0))=mean(diffarr)    ;eliminates pixels with errors >5% in IS test.
     diffarr(*,0:1)=mean(diffarr(*,3))               ;eliminates rows 0 & 1, too little signal.

     Dev(modno,naz,nt1,nt2)=stddev(diffarr,/double)

      ;capture the low state...
      if Dev(modno,naz,nt1,nt2) lt lowdev then begin
       lowdev=Dev(modno,naz,nt1,nt2)
       lowdiffarr=diffarr
       lowdn=dn
      endif


     nt2=nt2+1
    endfor	;T2

    nt1=nt1+1
   endfor	;T1

    minarr(naz)=min(Dev(modno,naz,*,*))   ;minimum DN/sec at each Azimuth
    print,sstr(naz+1),'/',sstr(nazs),'  Az = ',sstr(azi),'  Dev = ',sstr(min(Dev(modno,naz,*,*)))
     wait,0.1

   azi=azi+azstep
   if azi ge 360 then azi=azi-360
 endfor	;az

if pass ge 2 then goto,next_model	;finish after 2nd pass...


 ; create finer grid for 2nd pass...
 pass=pass+1
 index=where(dev eq min(dev))
 result=array_indices(dev,index(0))

 ;test to see if result is at the edge of the test range (i.e. not a contained minimum)
 upper_envelope=[nazs-1,nt1-1,nt2-1]
 upper_end=where(result(1:3) eq upper_envelope)
 lower_end=where(result(1:3) eq 0)
 if lower_end(0) ne -1 or upper_end(0) ne -1 then begin
  print,'***The result is at the end of the range...please re-define range.***'
  goto,result
 endif

 az1=azarr(result(1)-1)-az
 az2=azarr(result(1)+1)-az
 if az2 lt az1 then az2=az2+360
 azstep=azstep/factor1
 T1a=t1arr(result(2)-1)
 T1b=t1arr(result(2)+1)
 T1step=T1step/factor2
 T2a=t2arr(result(3)-1)
 T2b=t2arr(result(3)+1)
 T2step=T2step/factor2

 ;print out info from first pass...
 index=where(dev eq min(dev))
 result=array_indices(dev,index(0))
 print,'Indices (model,AZ,T1,T2)= ',result
 print,'Upper_indices (AZ,T1,T2)=                    ',upper_envelope
 print,'model= ',model(result(0))
 az_result=azarr(result(1))
 print,'Azimuth= ',sstr(az_result)
 t1=t1arr(result(2))
 print,'T1, tip= ',sstr(t1)
 t2=t2arr(result(3))
 print,'T2, tip= ',sstr(t2)
 ;ewtip=-(t1*cos((90-sza)*pi/180)*cos(esa*PI/180)+t2*cos((90-esa)*pi/180))
 ewtip=-(t1*cos(esa)+t2*sin(esa))
 print,'E/W tip= ',ewtip
 print,'Altitude (km)= ',alt
 print,'Estimated Az.= ',d_value(h,101)
 print,'Min Dev= ',lowdev,'  DN/sec, 1 sigma => ',100.*lowdev/mean(lowDN(*,*,ch1:ch2)),' %'

 goto,Pass


next_model:
endfor  ;model

;________________________________________
; Results...
      elapsed(4)=elapsed(4)+(systime(1)-oldt)/60.
      oldt=systime(1)


result:
index=where(dev eq min(dev))
result=array_indices(dev,index(0))

 ;test to see if result is at the edge of the test range (i.e. not a contained minimum)
 upper_envelope=[nazs-1,nt1-1,nt2-1]
 upper_end=where(result(1:3) eq upper_envelope)
 lower_end=where(result(1:3) eq 0)
 if lower_end(0) ne -1 or upper_end(0) ne -1 then begin
  print,'***The result is at the end of the range...please re-define range.***'
 endif

 print,'Indices (model,AZ,T1,T2)= ',result
 print,'Upper_indices (AZ,T1,T2)=                    ',upper_envelope
print,'model= ',model(result(0))
az_result=azarr(result(1))
print,'Azimuth= ',sstr(az_result)
t1=t1arr(result(2))
print,'T1, tip= ',sstr(t1)
t2=t2arr(result(3))
print,'T2, tip= ',sstr(t2)
;ewtip=-(t1*cos((90-sza)*pi/180)*cos(esa*PI/180)+t2*cos((90-esa)*pi/180))
ewtip=-(t1*cos(esa)+t2*sin(esa))
print,'E/W tip= ',ewtip
print,'Altitude (km)= ',alt
print,'Estimated Az.= ',d_value(h,101)
print,'Min Dev= ',lowdev,'  DN/sec, 1 sigma => ',100.*lowdev/mean(lowDN(*,*,ch1:ch2)),' %'

;Plot the final answer...

 SA_incident_intensities,model_name,alt,az,t1, $
   t2,intensity,dir


 DN=dnpersecV8(intensity,temps,wavl)	;DN(col,row,chain)

if ch1 eq 0 then data=(p(6:11,*)/(etime/1000.)) else data=(p(12:17,*)/(etime/1000.)) ;verticle
maxp=max([max(data),max(LowDN(*,*,ch1))])

window,0
Surface,LowDN(*,*,ch1),/save,zrange=[0,maxp],charsize=1.5,ztickformat='(i)',xtitle='Column',ytitle='Row',ztitle='DN/sec'
surface,data,/noerase,color=5000,zrange=[0,maxp],charsize=1.5,ztickformat='(i)'
xyouts,.5,.9,color+' Vert. Data (red) and Model (white), Az = '+sstr(az_result),alignment=0.5,/normal

if ch1 eq 0 then data2=(p(0:5,*)/(etime/1000.)) else data2=(p(18:23,*)/(etime/1000.)) ;horizontal
maxp=max([max(data2),max(LowDN(*,*,ch2))])

window,1
Surface,LowDN(*,*,ch2),/save,zrange=[0,maxp],charsize=1.5,ztickformat='(i)',xtitle='Column',ytitle='Row',ztitle='DN/sec'
surface,data2,/noerase,color=5000,zrange=[0,maxp],charsize=1.5,ztickformat='(i)'
xyouts,.5,.9,color+' Horiz. Data (red) and Model (white), Az = '+sstr(az_result),alignment=0.5,/normal

window,2
show3,lowdiffarr
xyouts,.1,.1,'Difference Array',alignment=0.5,/normal

window,3
plot,azarr,minarr,yrange=[0.9*min(minarr),1.4*min(minarr)],xrange=[fix(min(azarr)-1.5),fix(max(azarr)+2.5)]
xyouts,0.5,0.9,'Minimum DN/sec. at each Azimuth',alignment=0.5,/normal
xyouts,0.5,0.87,systime(0),alignment=0.5,/normal

exit:
total_time=(systime(1)-q)/60
print,'Processing Time = ',total_time,' minutes'
      elapsed(5)=elapsed(5)+(systime(1)-oldt)/60.
      oldt=systime(1)
print,"%time= Setup, Optimise1, SA_incident_intensities, DN/sec, Optmise2, Printing results"
print,100*elapsed/total_time,format='(6f21.2,"   %")'


stop
print,'Optimization Complete!'
end