c*** program to compute map from spherical harmonic coefficients c map value is written out in 5 degree increments c out.map contains the 5x5 degree map starting with as upper left c corner 87.5N,2.5E and ending with 87.5S 357.5E as lower right c corner. c out.map0 contains the 5x5 degree map starting with as upper left c corner 87.5N,177.5W and ending with 87.5S 177.5E as lower right c corner. c c This routine works with an input file that is formatted like c the bestcoff files. c-------FORMAT SPECIFICATIONS FOR BESTCOFF FILES----------------- c each mode has a set of lines: c line 1: n,l,and the number of s's (where mode: n_S_l, coefficient: c_s^t) c line 2: s for first s and 2*s+1 coefficients c line 3: 2*s+1 error values c c example: the first mode in bestcoffs is 0S2. c There are two s's available c_2^t and c_4^t. c c The bestcoff files only go up to s=6. c The c_s^t are in microHz and are ordered (for degree 2) c as c_2^0, Re(c_2^1), Im(c_2^1), Re(c_2^2), im(c_2^2) c and likewise for degree 4 c from c_4^0 to Im(c_4^4). c Using these cst's and this fortran code, the output maps c are in microHz. c c Change subroutine readmod to accomodate any other model format. c c The spherical harmonics are fully normalized (Edmonds, A.R. 1960, c Angular Momentum in Quantum Mechanics, Princeton Univ. Press). c c lmax.le.50!! c implicit real*8(a-h,o-z) parameter(lmx=51) common/map/amap(72,36) common/vmap/cof(2,lmx*(lmx+1)/2),lmax dimension x(101),x1(101),ss(36,101),cc(36,101) data pi/3.14159265358979d0/ rad=180.d0/pi call readmod(55,dum) open(66,file='out.map') open(67,file='out.map0') c*** compute sines and cosines only once do 25 m=1,lmax do 25 j=1,36 phi=(j*5.-2.5)/rad cc(j,m)=dcos(m*phi) 25 ss(j,m)=dsin(m*phi) c*** zero out map do 30 i=1,36 do 30 j=1,72 30 amap(j,i)=0. c*** accumulate map values, do all longitudes for each latitude c*** use symmetry to do 4 quadrants at once lp=lmax+1 do 50 i=1,18 ih=37-i theta=(i*5.-2.5)/rad co=dcos(theta) si=dsin(theta) ind=1 isn=-1 do 40 lp1=1,lp l=lp1-1 ind=ind+l isn=-isn call gxfcn(l,co,si,x,x1) do 45 j=1,36 jh=73-j ksn=isn amap(j,i)=amap(j,i)+x(1)*cof(1,ind) amap(jh,i)=amap(jh,i)+x(1)*cof(1,ind) amap(j,ih)=amap(j,ih)+ksn*x(1)*cof(1,ind) amap(jh,ih)=amap(jh,ih)+ksn*x(1)*cof(1,ind) do 45 m=1,l ksn=-ksn dum1=2.d0*x(m+1)*cof(1,ind+m)*cc(j,m) dum2=2.d0*x(m+1)*cof(2,ind+m)*ss(j,m) amap(j,i)=amap(j,i)+(dum1-dum2) amap(jh,i)=amap(jh,i)+(dum1+dum2) amap(j,ih)=amap(j,ih)+ksn*(dum1-dum2) 45 amap(jh,ih)=amap(jh,ih)+ksn*(dum1+dum2) 40 continue 50 continue smin=100000. smax=-100000. do 70 i=1,36 do 70 j=1,72 if(amap(j,i).gt.smax)smax=amap(j,i) if(amap(j,i).lt.smin)smin=amap(j,i) 70 continue print*,' min/max in map: ',smin,smax c*** write out map do 60 i=1,36 write(67,910) + (amap(j,i),j=37,72),(amap(j,i),j=1,36) 60 write(66,910) (amap(j,i),j=1,72) 910 format(30(8e15.5/)) stop end subroutine gxfcn(l,c,s,x,x1) c c=cos(theta), s=sin(theta),theta=colatitude,l=angular order, c x(1) contains m=0, x(2) contains m=1, x(k+1) contains m=k c m=azimuthal order 0.le.m.le.l . x1=dx/dtheta c x(l,m,theta)*exp(imphi) are a set of orthonormal spherical harmonics. implicit real*8(a-h,o-z) parameter (lmx=41) dimension x(1),x1(1),iexp(lmx),scl(6) data pi4/12.56637061435916d0/ data scl/1.d0,1.d-60,1.d-120,1.d-180,1.d-240,1.d-300/ lp1=l+1 fl2p1=l+lp1 con=dsqrt((fl2p1)/pi4) if(l.eq.0) then x(1)=con x1(1)=0.d0 return end if if(l.ge.lmx) then print *,'array dimension in gxfcn need expanding' stop end if c*** handle small arguments if(dabs(s).lt.1.d-20) goto 300 cot=c/s c*** try m recursions ; first compute xll f=1.d0 do 100 i=1,l 100 f=f*(i+i-1.d0)/(i+i) x(lp1)=con*dsqrt(f)*(-s)**l C*** use scaling if xll underflows if(dabs(x(lp1)).lt.1.d-295) goto 200 x1(lp1)=x(lp1)*l*cot c*** use m recurrence starting from m=l do 110 i=1,l m=lp1-i mp1=m+1 f=dsqrt(i*(fl2p1-i)) x(m)=-(x1(mp1)+m*cot*x(mp1))/f 110 x1(m)=(m-1)*x(m)*cot+x(mp1)*f return c c*** use m recurrence starting from m=l with scaling 200 x(lp1)=(-1.d0)**l x1(lp1)=x(lp1)*l*cot iexp(lp1)=0 do 210 i=1,l m=lp1-i mp1=m+1 iexp(m)=iexp(mp1) f=dsqrt(i*(fl2p1-i)) x(m)=-(x1(mp1)+m*cot*x(mp1))/f x1(m)=(m-1)*x(m)*cot+x(mp1)*f 211 if(abs(x(m)).gt. 1.d60) then iexp(m)=iexp(m)+1 x(m)=x(m)*1.d-60 x1(m)=x1(m)*1.d-60 goto 211 end if 210 continue c*** now unexponentiate and calculate sum sum=x(1)*x(1) do 225 i=2,lp1 iex=iexp(1)-iexp(i)+1 if(iex.lt.7) then x(i)=x(i)*scl(iex) x1(i)=x1(i)*scl(iex) sum=sum+2.d0*x(i)*x(i) else x(i)=0.d0 x1(i)=0.d0 end if 225 continue sum=con/dsqrt(sum) do 230 i=1,lp1 x(i)=x(i)*sum 230 x1(i)=x1(i)*sum return c c*** handle very small arguments 300 do 310 i=1,lp1 x(i)=0.d0 310 x1(i)=0.d0 x(1)=con x1(2)=-0.5d0*con*dsqrt(dble(l*lp1)) return end subroutine readmod(io,dum2) c read in cst file c lmax must be .lt. 50 c dimension of cof is cof(2,(lmax+1)(lmax+2)/2) c this routine sets odd harmonic degree to zero. implicit real*8(a-h,o-z) parameter(lmx=51) common/vmap/cof(2,lmx*(lmx+1)/2),lmax dimension rcof(2*lmx+1,lmx),rerr(2*lmx+1,lmx) character*256 filnam logical log1 40 continue print *,'enter filename for sph. harm. coeff. (*.lm):' read(*,'(a)') filnam inquire(file=filnam,exist=log1) if(.not.log1)goto 40 print*,'enter desired n and l' read(*,*)iin,iil open(io,file=filnam) 10 continue read(io,*,end=20) nn,ll,jj do i=1,jj lp=4*i+1 read(io,*)ij,(rcof(j,i),j=1,lp) read(io,*)(rerr(j,i),j=1,lp) enddo if(nn.ne.iin.or.ll.ne.iil)goto 10 lmax=min(lp,lmx-1) istop = (lmax+1)*(lmax+2)/2 do ind=1,istop cof(1,ind)=0 cof(2,ind)=0 enddo do i=1,jj ll=2*i ist=ll*(ll+1)/2+1 jj=1 cof(1,ist)=rcof(jj,i) cof(2,ist)=0. do j=1,ll ist=ist+1 jj=jj+1 cof(1,ist)=rcof(jj,i) jj=jj+1 cof(2,ist)=rcof(jj,i) enddo enddo close(io) return 20 continue print*,'desired mode not found: n,l= ',iin,iil close(io) stop end