pro OptimizeV3,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) 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.
;5) Run SA simulator at each point and calculate RMS deviation (data to simulator)
;    at each pixel.
;6) Keep the points bounding the minimum (RMS deviation) as endpoints and repeat 3
;    with a tighter grid.
;7) 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)

AzRangeToSun=40. ; Restriction on azimuth angle to the sun when selecting data.
az1=-5.	 ; negative excursion limit in Azimuth
az2=+8.         ; positive excursion limit in Azimuth
azstep=1.0       ; Initial azimuth step size...final is initial devided by 'factor'
T1a=-30.         ; Limits on T1, T2, Step_size, etc as with Azimuth above...
T1b=+15.
T1step=5.
T2a=-30.
T2b=+15.
T2step=5.
factor=10	     ; factor by which grid is reduce for second pass




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 (Rads), 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 AzRangeToSun 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)

 ;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


     ; First call takes longer because it must read the model
     SA_incident_intensities,model_name,alt,azi,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)

      ;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


endfor  ;model

;________________________________________
; Results...

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
print,'LowDev= ',lowdev

;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

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

window,2
show3,lowdiffarr

exit:
stop
end