FUNCTION div_iso,Vxx,Vyy,seuil ; ;programme de calcul de la "divergence isotropique" ;INPUTS: Vxx,Vyy - velocity field ; seuil - threshold (minimum value of iso.divergence - optional) ; ;CALLING: diverg = div_iso(vx,vy [,thr]) ; ;written by Thierry Roudier ;I/O and small things modified by Michal, 24.02.2003 ; norm = 2. norm2= sqrt(8) norm3= 4. norm4= sqrt(32) norm5= sqrt(20) taille = size(Vxx) nx=taille(1) & ny=taille(2) diverg=fltarr(nx,ny) ; Constants ; angles teta1 = (!pi/4) teta2 = 26.565051 *(!pi)/180 teta3 = 63.434949*(!pi)/180 ; sines and cosines ct1 = cos(teta1) & st1 = sin(teta1) ct2 = cos(teta2) & st2 = sin(teta2) ct3 = cos(teta3) & st3 = sin(teta3) ;----------------------------------------------------- for j=2,ny-3 do begin ; LOOP over Vyy points for i=2,nx-3 do begin ; LOOP over Vxx points ; 8 premiers proches voisins 1-ere couronne v0=Vxx(i+1,j) v1=Vxx(i+1,j+1) * ct1 + Vyy(i+1,j+1) *st1 v2= Vyy(i,j+1) v3=-Vxx(i-1,j+1) * st1 + Vyy(i-1,j+1) *ct1 v4= Vxx(i-1,j) v5=Vxx(i-1,j-1) * ct1 + Vyy(i-1,j-1) *st1 v6=Vyy(i,j-1) v7=-Vxx(i+1,j-1) * st1 + Vyy(i+1,j-1) *ct1 diver1= (v0-v4)/norm & diver2= (v1-v5)/norm2 diver3= (v2-v6)/norm & diver4= (v3-v7)/norm2 ; (divergence classique = diver1+diver3) ; 8 premiers proches voisins 2-eme couronne v00= Vxx(i+2,j) v11= Vxx(i+2,j+2) * ct1 + Vyy(i+2,j+2) * st1 v22= Vyy(i,j+2) v33=-Vxx(i-2,j+2) * st1 + Vyy(i-2,j+2) * ct1 v44= Vxx(i-2,j) v55= Vxx(i-2,j-2) * ct1 + Vyy(i-2,j-2) * st1 v66= Vyy(i,j-2) v77=-Vxx(i+2,j-2) * st1 + Vyy(i+2,j-2) * ct1 diver5= (v00-v44)/norm3 & diver6= (v11-v55)/norm4 diver7= (v22-v66)/norm3 & diver8= (v33-v77)/norm4 ; 8 autres proches voisins 2-eme couronne v000= Vxx(i+2,j+1) * ct2 + Vyy(i+2,j+1) * st2 ; vx v111= Vxx(i+1,j+2) * ct3 + Vyy(i+1,j+2) * st3 ; vx v222= -Vxx(i-1,j+2) * st2 + Vyy(i-1,j+2) * ct2 ; vy v333= -Vxx(i-2,j+1) * st3 + Vyy(i-2,j+1) * ct3 ; vy v444= Vxx(i-2,j-1) * ct2 + Vyy(i-2,j-1) * st2 ; vx v555= Vxx(i-1,j-2) * ct3 + Vyy(i-1,j-2) * st3 ; vx v666= -Vxx(i+1,j-2) * st2 + Vyy(i+1,j-2) * ct2 ; vy v777= -Vxx(i+2,j-1) * st3 + Vyy(i+2,j-1) * ct3 ; vy diver9= (v000-v444)/norm5 & diver10= (v111-v555)/norm5 diver11= (v222-v666)/norm5 & diver12= (v333-v777)/norm5 ; somdiv1 = diver1 + diver2 + diver3 + diver4 somdiv2 = diver5 + diver6 + diver7 + diver8 somdiv3 = diver9 + diver10 + diver11 + diver12 somdiv = somdiv1 + somdiv2 + somdiv3 ;filling the result array diverg(i,j)=somdiv endfor ; end LOOP Vxx endfor ; end LOOP Vyy ;thresholding the result to detect centers of positive divergence if n_params() gt 2 then begin ma=diverg gt seuil diverg=diverg*ma endif ; RETURN,diverg end