pro OptimizeV6,alt,ModNum

alt=80
ALT=88
alt=100
alt=60.4481	;test case
q=systime(1)   ;how long does it take?

;CSee 26Feb07, revised 07May2007

;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...

model_dir='Y:\Titan_ULVS_models\Model_777'  ;directory containing models
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)

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
  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

     ;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]

;____________________
;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)

 az1=-10.	 ; negative excursion limit in Azimuth
 az2=+15.         ; positive excursion limit in Azimuth
 azstep=5.0       ; 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=5.
 T2a=-30.
 T2b=+20.
 T2step=5.
 factor=10.	     ; factor by which grid is reduce for second pass

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


 ;____________________
 ;Optimization loop....

 print,'Optimize...'
 wait,0.5

 pass=1

 Pass:	;Two passes; first with course grid second with smaller, finer grid

  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)


  azi=az+az1	               ;start with lowest azimuth in range (offset from Erichs estimate)

  for naz=0,nazs-1 do begin    ;for each azimuth


   azarr(naz)=azi

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

    nt2=0	;index for T2
    for T2i=T2a,T2b,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
     SA_incident_intensities,model_name,alt,azimuth,tip_away_from_sun, $
       orthogonal_tip,intensity,dir

     ;;print,'calculating DN/sec...'
     ;;wait,0.5
     DN=dnpersec(intensity,temps,wavl)	;DN(col,row,chain)




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

     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

     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 ne -1 or upper_end 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/factor
 T1a=t1arr(result(2)-1)
 T1b=t1arr(result(2)+1)
 T1step=T1step/factor
 T2a=t2arr(result(3)-1)
 T2b=t2arr(result(3)+1)
 T2step=T2step/factor

 ;print out info from first pass...
 index=where(dev eq min(dev))
 result=array_indices(dev,index(0))
 print,'Indices= ',result,'   (model,AZ,T1,T2)'
 print,'Upper_indices =     ',upper_envelope,'     (AZ,T1,T2)'
 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'

 goto,Pass


next_model:
endfor  ;model

;________________________________________
; Results...

Result:
index=where(dev eq min(dev))
result=array_indices(dev,index(0))
print,'Indices= ',result
print,'Upper_indices =     ',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'

;Plot the final answer...

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


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

data=(p(6:11,*)/(etime/1000.))
maxp=max([max(data),max(LowDN(*,*,0))])

window,0
Surface,LowDN(*,*,0),/save,zrange=[0,maxp],charsize=1.5
surface,data,/noerase,color=5000,zrange=[0,maxp],charsize=1.5
xyouts,.5,.9,'Blue Vert. Data (red) and Model (white)',alignment=0.5,/normal

data2=(p(0:5,*)/(etime/1000.))
maxp=max([max(data2),max(LowDN(*,*,1))])

window,1
Surface,LowDN(*,*,1),/save,zrange=[0,maxp],charsize=1.5
surface,data2,/noerase,color=5000,zrange=[0,maxp],charsize=1.5
xyouts,.5,.9,'Blue Horiz. Data (red) and Model (white)',alignment=0.5,/normal

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

exit:
print,(systime(1)-q)/60,' minutes'
stop
end