;===============================================================
;pattern_main.pro
;Purpose: find the movement of smoke plume
;written by Jian Zeng
;2006-09-17
;2006-12-14       find the smoke movement only use two images
;===============================================================
  pro read_binary,outfile,data
        openr,1,outfile
        readu,1,data
        close,1
  end
;--------------------------------------------------------------------------
  pro find_move,target,search,glat,glon,outidx,$
                target_xs=target_xs,target_ys=target_ys,mxcor=mxcor


    pct=0.2

    index=where(target gt 0.,nin)
    outidx=[-1,-1]
    mxcor=-999.0

    if nin gt target_xs*target_ys*pct then begin

      step=1
      nx_sub=target_xs/step+1
      ny_sub=target_ys/step+1
      r=fltarr(nx_sub,ny_sub)-999.0


      for m=0,nx_sub-1,step do begin
        for n=0,ny_sub-1,step do begin
            i1=m
            i2=m+target_xs
            j1=n
            j2=n+target_ys

            sub=search(i1:i2,j1:j2)
            lat=glat(i1:i2,j1:j2)
            lon=glon(i1:i2,j1:j2)

            index2=where(lat gt -90 and lon gt -180. and sub gt 0.$
	                 and target gt 0.,nin2)
            if nin2 gt target_xs*target_ys*pct then begin
              
	      target_mean=mean(target(index2))
              target_std=total((target(index2)-target_mean)^2)
              sub_mean=mean(sub(index2))
              sub_std=total((sub(index2)-sub_mean)^2)
	      
              up=total((sub(index2)-sub_mean)*(target(index2)-target_mean))
	      if sqrt(target_std*sub_std) gt 0. then $
	        r(m,n)=up/sqrt(target_std*sub_std)

	    endif
	endfor
      endfor

      mxcor=max(r,location)
      index=where(abs(mxcor-r) lt 1.e-2,cp)
      if mxcor gt 0. then begin
        outidx=array_indices(r,location)
      endif

    endif ;if nin gt target_xs*target_ys*pct
    
end
;------------------------------------------------------------------
pro create_newimg,nx_g,ny_g,num,move,target_xs,target_ys,aod,newaod

 newaod=fltarr(nx_g,ny_g)-9999.0
 for inum=0,num-1 do begin
   ii0=move(inum,0)
   jj0=move(inum,1)
   if ii0 gt 0 and jj0 gt 0 then begin
     i1=(ii0-target_xs) > 0
     i2=(ii0+target_xs) < (nx_g-1)
     j1=(jj0-target_ys) > 0
     j2=(jj0+target_ys) < (ny_g-1)
     newaod(i1:i2,j1:j2)=aod(i1:i2,j1:j2)
   endif
 endfor

end
;-----------------------------------------------------------------------------
pro find_smokegrid,yr,month,nday,time_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
                   target_xs,target_ys,temp,nsmoke,data_suffix
    
    date_last=yr+month+nday+time_last
    infile_aod=file_search(dir0+'*'+date_last+data_suffix+'*',count=nfl_aod)
;    infile_aod=file_search(dir0+'/*'+date_last+data_suffix+'*',count=nfl_aod)
    print,'in find_smokegrid'+dir0+'/*'+date_last+data_suffix+'*'
    maxnum=nx*ny
    temp=lonarr(maxnum,2)
    nplumes=0
    plumes=lonarr(maxnum,2)
    nsmoke=0L
    npix=target_xs*target_ys*0.  ;0.3 changed on Jan. 15, 2008

    if nfl_aod eq 1 then begin
      saod=fltarr(nx_g,ny_g)
      read_binary,infile_aod,saod
      
      for inx=0,nx-1 do begin
      for iny=0,ny-1 do begin
        mini=origx(inx)
        minj=origy(iny)
        i1=(mini-target_xs/2) > 0
        i2=(mini+target_xs/2) < (nx_g-1)
        j1=(minj-target_ys/2) > 0
        j2=(minj+target_ys/2) < (ny_g-1)
        target=saod(i1:i2,j1:j2)
        ind=where(target gt 0.,cp)
        if cp gt npix then begin
          temp(nsmoke,0)=mini
          temp(nsmoke,1)=minj
          nsmoke=nsmoke+1
        endif
      endfor
      endfor
      print,'Smoke Grids (#, max, min) at ',date_last,' = ',$
             nsmoke,max(saod),min(saod(where(saod gt -9)))
    endif
    
end
;================================================================
;main program
;================================================================
pro pattern_main,dir,dir0,yr,jday,time,nx_g,ny_g,newaod,glat,glon,data_suffix
              
  target_xs=32
  target_ys=32
  nx=fix(nx_g/target_xs)-2
  ny=fix(ny_g/target_ys)-2
    
  file=file_search(dir+yr+jday+'*all.aod.gz',count=nfl0) 
  print,'file in pattern dir'+dir+yr+jday+'*all*'
  
  if nfl0 eq 0 then $
    print, 'Pattern_main.pro: No AOD file found on the day '+jday+' of the year '+yr $
  else begin
  
    convert_julianday,yr,jday,mon,day
    month=string(mon,format='(i2.2)')
    nday=string(day,format='(i2.2)')
    date0=yr+month+nday+time
    print, 'starting from ',date0
    pos=strpos(file(0),'_i18_US.all.aod.gz')
    times=strmid(file,pos-6,4)
    part=where(times eq time,nfl)

    if nfl ne 1 then print,'no AOD file found at '+time+'Z' $
    else if part(0) eq 0 then $
      print,'no previous AOD file can be found' $
      
    else begin
      infile=file(part)
      infile_last=file(part-1)
 
      origx=indgen(nx)*target_xs+target_xs*3/2
      origy=indgen(ny)*target_ys+target_ys*3/2
      maxnum=nx*ny
      
      move=lonarr(maxnum,2)-9999
      corr=fltarr(maxnum)-9999.0
  
      time_last=strmid(infile_last,pos-6,4)
      date_last=yr+month+nday+time_last

      print,'in pattern dir0'+dir0
      find_smokegrid,yr,month,nday,time_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
                     target_xs,target_ys,temp,nsmoke,data_suffix

      if nsmoke eq 0L and part(0) gt 1 then begin
        
        print, 'no smoke grid found 30 minutes ago, search one hour ago...'
       
        infile_last=file(part-2)
        time_last=strmid(infile_last,pos-6,4)
        date_last=yr+month+nday+time_last
        find_smokegrid,yr,month,nday,time_last,dir0,nx,ny,origx,origy,nx_g,ny_g,$
                       target_xs,target_ys,temp,nsmoke,data_suffix

      endif  ;if nsmoke eq 0L and part(0) ge 1 

      if nsmoke eq 0L then begin
        print,'no smoke grid found within past 60 minutes' 
      endif else begin
        
        print,'find pattern at -- ',date_last
        print,'in pattern infile '+infile
        print,'in pattern         aod=read_GASP('+infile_last
        aod0=read_GASP(infile,nx_g=nx_g,ny_g=ny_g)

        aod=read_GASP(infile_last,nx_g=nx_g,ny_g=ny_g)

        ind=where(glat gt -90. and glon gt -180.)

        for inum=0,nsmoke-1L do begin
          mini=temp(inum,0)
          minj=temp(inum,1)
          if mini ge origx(0) and minj ge origy(0) and $
	     mini le origx(nx-1) and minj le origy(ny-1) then begin
            i1=(mini-target_xs/2) > 0
            i2=(mini+target_xs/2) < (nx_g-1)
            j1=(minj-target_ys/2) > 0
            j2=(minj+target_ys/2) < (ny_g-1)

            target=aod0(i1:i2,j1:j2)

	    i1=(mini-target_xs) > 0
            i2=(mini+target_xs) < (nx_g-1)
            j1=(minj-target_ys) > 0
            j2=(minj+target_ys) < (ny_g-1)
            search=aod(i1:i2,j1:j2)

            find_move,target,search,glat,glon,outidx,$
                      target_xs=target_xs,target_ys=target_ys,$
                      mxcor=mxcor

            if outidx(0) ne -1 then begin
              move(inum,0)=(outidx(0)-target_xs/2+mini) $
	          > (origx(0)-target_xs/2) < (origx(nx-1)+target_xs/2)
              move(inum,1)=(outidx(1)-target_ys/2+minj)$
	          > (origy(0)-target_ys/2) < (origy(ny-1)+target_ys/2)
            endif else begin
              move(inum,0)=mini
              move(inum,1)=minj
	    endelse

            corr(inum)=mxcor

          endif else begin ;if mini gt 0 and minj gt 0
            move(inum,0)=mini
            move(inum,1)=minj

          endelse
        endfor  ;for inum=0,num-1L do begin

        create_newimg,nx_g,ny_g,nsmoke,move,target_xs,target_ys,aod0,newaod
       
      endelse ;if nsmoke eq 0L 
      
    endelse ;if nfl ne 1 or part(0) eq 0 
    
  endelse ;if nfl0 eq 0 

end