pro diverg_lap_petits,seuil ;adapte aux donnees de Lapalma le 30 Aout 2001 ;programme de calcul de la divergence ; pour les donnees 15min (pas 15 min) ; boite 12 pixels (1.5 arcsec) ; donnees sequence 3 heure du Pic du midi ; verifie le 7 Septembre 1999 fenetre=5 ; fenetre de 5 pixels pour le calcul des vitesses moyennes. norm =2 * (0.125 * 736. * 12.) norm2=sqrt(8)*.125 * 736. * 12. norm3=4*(0.125 * 736. * 12.) norm4=sqrt(32)*.125 * 736. * 12. norm5=sqrt(20)*.125 * 736. * 12. nt=99 ;for ii=0,nt-1 do begin for ii=0,32 do begin if( ii lt 9 )then begin fich=strcompress(string(format='("vipx",i1,".fits")',ii+1)) print,fich vxx=readfits(fich) fich2=strcompress(string(format='("vipy",i1,".fits")',ii+1)) print,fich2 vyy=readfits(fich2) endif else begin fich=strcompress(string(format='("vipx",i2,".fits")',ii+1)) print,fich vxx=readfits(fich) fich2=strcompress(string(format='("vipy",i2,".fits")',ii+1)) print,fich2 vyy=readfits(fich2) endelse ; (728 * 0.14arcsec/20sec)=5.096 Km/sec ;(736 * 0.125 arcsec/21sec)=4.381 Km/sec Vxx= Vxx * 4.381 ; normalisation en Km/sec apres la LCT Vyy= Vyy * 4.381 ; pour les donnees de Lapalma 9h ;Vxx= Vxx * 5.096 ; normalisation en Km/sec apres la LCT ;Vyy= Vyy * 5.096 ; pour les donnees du Pic du Midi 3h taille = size(Vxx) nx=taille(1) & ny=taille(2) print,'nx,ny',nx,ny diverg=fltarr(nx,ny) divrai=fltarr(nx,ny) ;rotation=fltarr(nx,ny) print,'min et max des vitesses vx',min(Vxx),max(Vxx) print,'min et max des vitesses vy',min(Vyy),max(Vyy) ; angle pour la 1ere couronne teta1= (!pi/4) for j=2,ny-3 do begin for i=2,nx-3 do begin ; 8 premiers proches voisins 1ere couronne v0=Vxx(i+1,j) v1=Vxx(i+1,j+1) * cos(teta1) + Vyy(i+1,j+1) *sin(teta1) v2= Vyy(i,j+1) v3=-Vxx(i-1,j+1) * sin(teta1) + Vyy(i-1,j+1) *cos(teta1) v4= Vxx(i-1,j) v5=Vxx(i-1,j-1) * cos(teta1) + Vyy(i-1,j-1) *sin(teta1) v6=Vyy(i,j-1) v7=-Vxx(i+1,j-1) * sin(teta1) + Vyy(i+1,j-1) *cos(teta1) diver1= (v0-v4)/norm & diver2= (v1-v5)/norm2 diver3= (v2-v6)/norm & diver4= (v3-v7)/norm2 ; print,'divergence classique =====', diver1+diver3 ; 8 premiers proches voisins 2eme couronne v00= Vxx(i+2,j) v11= Vxx(i+2,j+2) * cos(teta1) + Vyy(i+2,j+2) * sin(teta1) v22= Vyy(i,j+2) v33=-Vxx(i-2,j+2) * sin(teta1) + Vyy(i-2,j+2) * cos(teta1) v44= Vxx(i-2,j) v55= Vxx(i-2,j-2) * cos(teta1) + Vyy(i-2,j-2) * sin(teta1) v66= Vyy(i,j-2) v77=-Vxx(i+2,j-2) * sin(teta1) + Vyy(i+2,j-2) * cos(teta1) diver5= (v00-v44)/norm3 & diver6= (v11-v55)/norm4 diver7= (v22-v66)/norm3 & diver8= (v33-v77)/norm4 ; 8 autres proches voisins 2eme couronne teta2=26.565051 *(!pi)/180 teta3= 63.434949*(!pi)/180 v000= Vxx(i+2,j+1) * cos(teta2) + Vyy(i+2,j+1) * sin(teta2) ; vx v111= Vxx(i+1,j+2) * cos(teta3) + Vyy(i+1,j+2) * sin(teta3) ; vx v222= -Vxx(i-1,j+2) * sin(teta2) + Vyy(i-1,j+2) * cos(teta2) ; vy v333= -Vxx(i-2,j+1) * sin(teta3) + Vyy(i-2,j+1) * cos(teta3) ; vy v444= Vxx(i-2,j-1) * cos(teta2) + Vyy(i-2,j-1) * sin(teta2) ; vx v555= Vxx(i-1,j-2) * cos(teta3) + Vyy(i-1,j-2) * sin(teta3) ; vx v666= -Vxx(i+1,j-2) * sin(teta2) + Vyy(i+1,j-2) * cos(teta2) ; vy v777= -Vxx(i+2,j-1) * sin(teta3) + Vyy(i+2,j-1) * cos(teta3) ; 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 ; divergence pour reperer les explosions isotropes ; if (somdiv gt seuil) then diverg(i,j)=somdiv ; ; divergences classiques divrai(i,j)=diver1+diver3 endfor endfor print,min(diverg),max(diverg) print,'min et max des divergences en sec**-1',min(diverg),max(diverg) ; divergences vraies window,0 bob=rebin(divrai,420,420) tvscl,bob window,1 ; ;divergence pour reperer les explosions isotropes ; b=rebin(diverg,420,420) val=size(b) print,'taille de b',val(1),val(2) tvscl,b ; bib=b ; modif le 13 juillet 99 pour utilisation de cgrav_surf_final_blobs if( ii lt 9 )then begin fich2=strcompress(string(format='("div15_lap_pet",i1)',ii+1));ww version 18 12 1999 save,bib,file=fich2 endif else begin fich2=strcompress(string(format='("div15_lap_pet",i2)',ii+1));wa version 18 12 1999 save,bib,file=fich2 endelse if( ii lt 9 )then begin fich3=strcompress(string(format='("div15_lapvr_pet",i1)',ii+1));vr version 18 12 1999 save,bob,file=fich3 endif else begin fich3=strcompress(string(format='("div15_lapvr_pet",i2)',ii+1));vr version 18 12 1999 save,bob,file=fich3 endelse ; endfor ; end