source: trunk/src/interpolswath.pro @ 583

;
;+
;
; .. _interpolswath.pro:
;
; =================
; interpolswath.pro
; =================
;
; DESCRIPTION
; ===========
;
; methode d'affectation des tb amsu sur grille de resolution
; fine pour eviter les pbs de trous entre les pixels
;
; tb du canal - dimension fova (30 points pour amsua)
; latlu
; lonlu memes dimensions que tb
;
; principe: on lit la fauchee, puis on calcule fova, distances de puis
; le centre de chaque pixel (en principe symetrique)
; on cree grid, de dimension ngrid pixels, depuis le centre
; jusqu'au bout du dernier pixel de chaque cote, avec une
; resolution fixee au diametre du pixel central ou la moitie (resol
; vaut 1 ou 2)
;
; on interpole les longitudes et latitudes en correspondance
; on affecte a chaque pixel nouveau la Tb originelle la plus proche
; en sortie on a une fauchee synthetique amsua reguliere
; TODO
; =====
;
; lelod 20120530
; faire evoluer le prgm pour l'inetrpolation des donnees AMSUB
;
; regler le pb des longitudes autour de 180 degres
; l'interpolation ne les gere pas et cree des points dans la zone AMMA...
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20120417
; 
;   * usage of key_performance
;
; - fplod 20111130T120358Z cratos (Linux)
;
;   * minimal rest header
;
;-
pro interpolswath, tb, latlu, lonlu, masklu, resol, nbgrid, tbint, latgrid, longrid, mask

@common
s = size(SCOPE_TRACEBACK(/STRUCTURE),/DIMENSION)
routine = (SCOPE_TRACEBACK(/STRUCTURE))[s-1].Routine
;
if key_performance EQ 1 THEN BEGIN
    time1 = SYSTIME(1)
ENDIF
;
pixelsize,pixatot,pixbtot,alongatot,alongbtot
na=max(size(pixatot))
nnadir=na/2
pixa=pixatot[nnadir:(na-1)]
fova=fltarr(na)
fova[nnadir]=pixa[0]/2.
fova[nnadir-1]=-pixa[0]/2.
rterre=6400.
deg=!pi/180.
lat=latlu*deg
lon=lonlu*deg
mnadir=nnadir-1
ifov=indgen(na)
for i=nnadir+1,na-1 do begin
   dist=rterre*acos(sin(lat(i-1))*sin(lat(i)) + cos(lat(i-1))*cos(lat(i))*cos(lon(i)-lon(i-1)))
   fova[i]=fova[i-1]+dist
endfor
for i=nnadir-2,0,-1 do begin
   dist=rterre*acos(sin(lat(i))*sin(lat(i+1)) + cos(lat(i))*cos(lat(i+1))*cos(lon(i)-lon(i+1)))
;print,i,dist
   fova[i]=fova[i+1]-dist
endfor
;print,'fova'
;print,fova
;print,'pixatot'
;print,pixatot
;window,1
 ;ploterr,ifov,fova,pixatot

distmax=-min(fova)+max(fova)+max(pixa)
dmax=distmax/2.
delta= fix(pixa[0]/resol)
;print,'pas swath',delta
ndemigrid=fix(dmax/delta)+1
if (dmax/delta-fix(dmax/delta) le 0.5) then ndemigrid=ndemigrid-1
nbgrid=ndemigrid*2
;print,'nb pts',nbgrid, ' dmax',dmax

grid2=indgen(ndemigrid)*delta+fova[nnadir]
grid1=fltarr(ndemigrid)
for i=0,ndemigrid-1 do begin
   grid1[i]=-grid2[ndemigrid-1-i]
endfor
;alternative equivalente)
;for i=ndemigrid-1,0,-1 do grid1[i]=-(ndemigrid-i-1)*delta-fova[nnadir]
grid=[grid1,grid2]
;print,grid
;print,tb
tbint=fltarr(nbgrid)
mask=intarr(nbgrid)
for i=0,nbgrid-1 do begin
   ind=where(abs(grid[i]-fova) eq min(abs(grid[i]-fova)),nii)
      if nii eq 1 then begin
         tbint[i]=tb[ind]
         mask[i]=masklu[ind]
      endif else begin
         print,'pas de minimum?',nii
      endelse
   endfor

;for i=0,na-1 do begin
 ;  ind=where(abs(grid-fova[i]) le pixatot[i]/2,nii)
  ; print,i,nii
  ; if nii ne 0 then begin
   ;   tbint[ind]=tb[i]
  ;    mask[ind]=masklu[i]
  ; endif
  ; if nii eq 0 then begin
  ;     
;endfor


latgrid=interpol(latlu,fova,grid)
longrid=interpol(lonlu,fova,grid)

;print,tbint
IF key_performance EQ 1 THEN BEGIN
   msg = report(['ppp : ' + routine + ' : durée totale ' $
               + string(SYSTIME(1)-time1,format='(F12.6)')])
ENDIF
;
end