pro Optimize,alt,ModNum

;CSee 26Feb07
;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) Make a list of the parameters (specifically altitude, but also temperatures,
;     azimuths, and E/W tip) from the available Titan descent data.
;2) Define initial point using Az from header, Ti=T2=0
;3) Create a grid of values using a range of +/- 15 deg and 5 points.
;4) Keep the two bounding the minimum (RMS deviation) as range endpoints and repeat 3
;5) Interpolate between two points bounding the min and compare to the min in 4,
;    use the least of these two. (10 calls to model)
;6) Repeat in this order: T1,T2,Az (30 calls to model)
;7) Repeat 2 thru 6 with new start point (60 calls total)
;8) Save Model, Az, T1, T2, Min_RMS_dev
;9) Repeat for next model.


alt=60.4481	;test case

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



common saved,model_name_old,alt_old,tables_read,wavl,solar_zen_angle, $
  theta,phi,idnc,zeniths_blue_vertical,azimuths_blue_vertical, $
  zeniths_blue_horizontal,azimuths_blue_horizontal

; 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, angle of the sun relative to east (114-90)




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

 ;find which Descent data to use...
 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 40 then goto,nextfile
  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

;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

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)

nazs=14
nt1s=10
nt2s=10
 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)

 ;____________________
 ;Optimization loop....

print,'Optimize...'
wait,0.5

 naz=0	;index for Az
 for azi=az-5,az+8,1 do begin   ;5 degrees each side of Erichs estimate


 nt1=0	;index for T1
 for T1i=T1-30,T1+15,5 do begin
 t1arr(nt1)=t1i

 nt2=0	;index for T2
 for T2i=T1-30,T2+15,5 do begin
 t2arr(nt2)=t2i

 tip_away_from_sun=T1i
 orthogonal_tip=T2i
 if azi lt 360 then azimuth=azi else azimuth=azi-360
 azarr(naz)=azimuth

 ;;probe_tip=0.0

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


 ; First call takes longer because it must read the model
 SA_incident_intensities,model_name,alt,azimuth,tip_away_from_sun, $
   orthogonal_tip,intensity,dir


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


 ; 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

 diffarr=[DN(*,*,0)-(p(6:11,*)/(etime/1000.)),DN(*,*,1)-(p(0:5,*)/(etime/1000.))]
 pixelmask=[pixelmaskbv,pixelmaskbh]
 diffarr(where(pixelmask eq 0))=mean(diffarr)

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

 nt2=nt2+1
 endfor	;T2

 nt1=nt1+1
 endfor	;T1

   print,'Az = ',sstr(azimuth),'  Dev = ',sstr(min(Dev(modno,naz,*,*)))
   wait,0.1
 naz=naz+1
endfor	;az


endfor  ;model

index=where(dev eq min(dev))
result=array_indices(dev,index(0))
print,result
print,'model= ',model(result(0))
az=azarr(result(1))
print,'Azimuth= ',sstr(az)
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



exit:
stop
end