;-------------------------------------------------------------------------------- ; This program takes the point-to-point Born kernels of Aaron C. Birch ; and change their ; format so that they can be used with my inversion code inversionRLS.pro ; It computes the average over an annulus IN A WAY COMPATIBLE WITH THE AVERAGING ; SCHEME OF CORRELATION.PRO ; and it divides the value by the sound speed at each depth ; ASSUMPTIONS AND SIMPLIFICATIONS: ; -The kernels are symmetrical around their center ; -We make a location error (when we shift) of less than half a pixel ; -We make a location error (when we rebin) ; -The sound-speed profile c(z) is assumed to be "model.fits" ; extracted from the file /home/aaronb/projects/kc/model.mat ; -Aaron's kernels are assumed to be in the format 800*200*100 ; NB: THE AZIMUTHAL AVERAGING IS MADE WITH BILINEAR INTERPOLATION ;--------------------------------------------------------------------------------- PRO BornAvgKern nkern = 11 ; we compute 11 azimuthally averaged kernels kernname = STRARR(nkern,5) ; each one is an average over 5 kernels (1 pixel-wide) savename = STRARR(nkern) abscisse = STRARR(nkern) ordonnee = STRARR(nkern) distances = FLTARR(nkern,5) FOR i=0l,nkern-1 DO BEGIN FOR j = 0l,4 DO kernname[i,j]='/scr109/couvidat/INVERSIONS/BORN/Bornavgkern'+STRTRIM(STRING(i+1),1)+STRTRIM(STRING(j+1),1)+'_new.fits' abscisse[i] ='/scr109/couvidat/INVERSIONS/BORN/xBornavg'+STRTRIM(STRING(i+1),1)+STRTRIM(STRING(1),1)+'_new.fits' ; same x for 5 annuli ordonnee[i] ='/scr109/couvidat/INVERSIONS/BORN/oBornavg'+STRTRIM(STRING(i+1),1)+STRTRIM(STRING(1),1)+'_new.fits' ; same y for 5 annuli savename[i] ='/scr109/couvidat/INVERSIONS/BORN/Bornavgk_'+STRTRIM(STRING(i+1),1)+'.bin' ENDFOR z = MRDFITS('/scr109/couvidat/INVERSIONS/BORN/zBornavg_new.fits') ; DISTANCE SOURCE-RECEIVER IN Mm distances[0,*] = 3.7+(8.7-3.7) /4.*FINDGEN(5) distances[1,*] = 6.2+(11.2-6.2) /4.*FINDGEN(5) distances[2,*] = 8.7+(14.5-8.7) /4.*FINDGEN(5) distances[3,*] = 14.5+(19.4-14.5)/4.*FINDGEN(5) distances[4,*] = 19.4+(29.3-19.4)/4.*FINDGEN(5) distances[5,*] = 26.0+(35.1-26.0)/4.*FINDGEN(5) distances[6,*] = 31.8+(41.7-31.8)/4.*FINDGEN(5) distances[7,*] = 38.4+(47.5-38.4)/4.*FINDGEN(5) distances[8,*] = 44.2+(54.1-44.2)/4.*FINDGEN(5) distances[9,*] = 50.8+(59.9-50.8)/4.*FINDGEN(5) distances[10,*] = 56.6+(66.7-56.6)/4.*FINDGEN(5) c = MRDFITS('/scr109/couvidat/INVERSIONS/BORN/model.fits') ; sound speed from the solar model used ; by Aaron to produce his kernels zlim = 300 ; I cut all the depths below !!!!WARNING: CHANGE INVERSIONRLS.PRO TOO IF YOU CHANGE THIS VALUE z = z[0:zlim-1] c = c[0:zlim-1] c = REVERSE(c) z = DOUBLE(REVERSE(z)) ; VERTICAL SCALE WE WANT FOR THE AVERAGED KERNELS z2 = [z[0],z[39],z[79],z[99],z[119],z[139],z[159],z[179],z[199],z[219],z[239],z[259],z[279]] samplingz = [z2[1:12]-z2[0:11],z[299]-z[279]] zelem = N_ELEMENTS(z2) ; NUMBER AND VALUES OF THE ANGLES WE USE FOR THE AZIMUTHAL AVERAGE nangles = 1000 angles = FINDGEN(nangles)/FLOAT(nangles-1)*360.0-180.0 ; we average each kernel over an annulus: ; BORN KERNELS MUST BE AVERAGED THE SAME WAY WE COMPUTE ; THE TRAVEL-TIME OVER AVERAGED ANNULI ;----------------------------------------------------------------------------------------------------- PRINT,'READING THE ABSCISSAE AND AVERAGING THE KERNELS' WINDOW,0,retain=2,xsize=600,ysize=600 !p.multi=0 FOR i=0,nkern-1 DO BEGIN y = MRDFITS(ordonnee[i]) dy = y[1]-y[0] ny = N_ELEMENTS(y) x = MRDFITS(abscisse[i]) dx = x[1]-x[0] nx = N_ELEMENTS(x) KernM = FLTARR(nx-1,nx-1,zelem) PRINT,'Resolutions: dx, dy, dy/dx,',dx,dy,dy/dx PRINT,'WARNING: IF dy/dx IS TOO DIFFERENT FROM 1 THEN THIS CODE HAS TO BE MODIFIED' FOR j=0,4 DO BEGIN PRINT,'ANNULUS INDEX',j Kern = MRDFITS(kernname[i,j]) Kern = Kern[0:zlim-1,*,*] Kern = REVERSE(Kern,1) distance = distances[i,j] IF(nx MOD 2 EQ 0) THEN x0 = (nx-1.)/2. - distance/2./dx ELSE x0 = nx/2. - distance/2./dx PRINT,'CENTER x0=',x0 ; NB: x0 IS NOT ROUNDED SO THAT ROT INTERPOLATES ; WE DIVIDE BY THE SOUND-SPEED ; BECAUSE AARON'S RAW KERNELS ARE MULTIPLIED BY c FOR k=0,zlim-1 DO Kern[k,*,*] = Kern[k,*,*]/c[k] Kern_z = FLTARR(zelem,nx,ny) ; WE INTERPOLATE THE KERNELS ON THE GRID IN DEPTH USED FOR THE ; INVERSION PROGRAM ;---------------------------------------------------------------------------------------------- PRINT,'INTEGRATION IN DEPTH' FOR iii=0,nx-1 DO BEGIN PRINT,iii FOR jjj=0,ny-1 DO BEGIN FOR zz=0l,zelem-1 DO BEGIN IF(zz NE zelem-1) THEN a = WHERE(z GE z2[zz] AND z LE z2[zz+1]) ELSE a = WHERE(z GE z2[zz]) IF (N_ELEMENTS(a) GT 1) THEN Kern_z[zz,iii,jjj]=INT_TABULATED(z[a],Kern[a,iii,jjj])/ABS(samplingz[zz]) ELSE PRINT,'PROBLEM !!!!',zz ; the discretization scheme is piecewise constant ENDFOR ENDFOR ENDFOR TVIM,Kern_z[0,*,*],/scale DELVAR,Kern ; to free space ; AZIMUTHAL AVERAGING (PLEASE: READ ALL THE ASSUMPTIONS) ;------------------------------------------------------------------- PRINT,'AVERAGING' kernel = FLTARR(nx-1,nx-1,zelem) ; we want an odd number of gridpoints to have a center FOR k=0,zelem-1 DO BEGIN ; we assume ny and nx are even numbers PRINT,k ; and we assume dy/dx=1.005 is actually =1, thus making a location error of 0.5% temp_re = FLTARR(nx-1,2*ny-1) temp_re[*,ny-1:2*ny-2] = REFORM(kern_z[k,1:nx-1,*]) temp_re[*,0:ny-1] = REVERSE(REFORM(kern_z[k,1:nx-1,*]),2) temp = FLTARR(nx-1,nx-1) FOR iii=0,nangles-1 DO temp = temp + ROT(temp_re,angles[iii],1,x0-1,ny-1,/INTERP) temp = temp /nangles;(dtheta*nangles) kernel[*,*,k] = temp TVIM,kernel[*,*,k],/scale; NB: this kernel is ignored in inversionRLS.pro ENDFOR ; SHIFT OF THE SOURCE POINT TO (0,0) Kernel = SHIFT(Kernel,-(nx/2-1),-(ny-1),0) ; AVERAGE OF 5 KERNELS KernM = kernel + KernM ENDFOR KernM = KernM/5.0 xdata = FINDGEN(nx-1)*dx xdata = xdata-xdata[nx/2-1] xdata = SHIFT(xdata,-(nx/2-1)) ydata = xdata ; WE ASSUME THAT dy=dx AND WE JUST SAVE dx SAVE,KernM,dx,z2,samplingz,xdata,ydata,file=savename[i] ENDFOR END