      program modeval 

      implicit real*8(a-h,o-z)
      dimension amap(72,36)
 
      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)



c  top of processing loop

        call readmod(7)
        write(6,*)'enter lmax for map'
        read(5,*)lmax2
	write(6,'(a)')'enter depth'
	read(5,*)rad
        rad=6371.d0-rad

	open (9,file='outmap')
      do 10 i=1,36
      do 10 j=1,72
       flat=92.5-i*5
       flong=-2.5+j*5
  10   call evmap(flat,flong,rad,amap(j,i),kopt,lmax2)
      if(kopt.ne.3) then
         if(kopt.eq.4) then
           anorm=-1.
         else
           print *,'enter average velocity'
           read(16,*) anorm
           if(kopt.eq.2) anorm=2.*anorm*anorm
         end if
          do 15 i=1,36
          do 15 j=1,72
   15       amap(j,i)=amap(j,i)/anorm
      end if
      do 20 i=1,36
  20     write(9,999)(amap(j,i),j=1,72)
  999 format(30(8e15.5/))
      close(9)
	write(6,*)'aspherical model slice stored in outmap'

      stop
      end

      subroutine evmap(flat,flong,radius,val,kopt,lmax2)
c  evaluates the perturbation in shear velocity, val,
c  at latitude flat(in deg), longitude flong(in deg) and radius (in km)
c  if kopt.eq.1 then perturbation is dv in km/s
c  if kopt.eq.2 then perturbation is d(v**2) or 2vdv in km/s
c  if kopt.eq.3 then perturbation is relative: dv/v
c  if kopt.eq.4 then perturbation is relative slowness: du/u
      implicit real*8(a-h,o-z)
      dimension x(51),x1(51),ss(51),cc(51),func(41)
      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)
      data pi/3.14159265358979d0/
      rad=180.d0/pi                                   
      do 5 i=1,nlay
        ilay=i
      if (iopt(i).ne.5) then
	rad1=radii(1,i)
	rad2=radii(2,i)
      else
        rad1=radii(3,i)
        rad2=radii(nfun(ilay)+2,i)
      endif
        if(radius.ge.rad1.and.radius.le.rad2) goto 10
    5   continue
      val=9999.
      return
   10 theta=(90.-flat)/rad
      phi=flong/rad
      do 25 m=1,lmax2
          cc(m)=dcos(m*phi)
   25     ss(m)=dsin(m*phi)        
      co=dcos(theta)
      si=dsin(theta)
      call radfunx(func,radius,ilay)
      kopt=ipow(ilay)
      val=0.
      ind=1
      ipar=iopt(ilay)
      do 40 lp1=1,lmax2+1
        l=lp1-1
        ind=ind+l
        call gxfcn(l,co,si,x,x1)
        if(ipar.eq.3) then
          rp=func(1)*cof(1,ind,ilay)+func(2)*cof(1,ind,ilay+1)
        else
          rp=0.
          do 41 k=1,nfun(ilay)
   41       rp=rp+func(k)*cof(2*k-1,ind,ilay)
        end if
        val=val+x(1)*rp
        do 45 m=1,l
        im=ind+m
        if(ipar.eq.3) then
          rp=func(1)*cof(1,im,ilay)+func(2)*cof(1,im,ilay+1)
          ap=func(1)*cof(2,im,ilay)+func(2)*cof(2,im,ilay+1)
        else
          rp=0.
          ap=0.
          do 44 k=1,nfun(ilay)
            rp=rp+func(k)*cof(2*k-1,im,ilay)
   44       ap=ap+func(k)*cof(2*k,im,ilay)
        end if
   45   val=val+2.d0*x(m+1)*(rp*cc(m)-ap*ss(m))
   40 continue
      return
      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 
c    harmonics.
      implicit real*8(a-h,o-z)
      dimension x(1),x1(1)
      data pi4/12.56637061435916d0/
      lp1=l+1
      fl2p1=l+lp1
      con=dsqrt((fl2p1)/pi4)
      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)=dsqrt(fl2p1*f/pi4)*(-s)**l
C*** use l recursions 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*** resort to l recursions ; first compute Pmm's
  200 x(1)=1.d0
      do 210 i=1,l
  210 x(i+1)=x(i)*s*(i+i-1.d0)
c*** perform l recurrence for each m value
      do 220 mm=1,l
      m=mm-1
      pmm=0.d0
      do 220 ll=mm,l
      pll=(c*(ll+ll-1.d0)*x(mm)-(ll+m-1.d0)*pmm)/(ll-m)
      pmm=x(mm)
  220 x(mm)=pll
c*** normalise and use m recurrence to compute derivative
      x(1)=x(1)*con
      x1(1)=0.d0
      if(l.gt.0) x1(1)=-x(2)*con
      do 230 mm=2,lp1
      m=mm-1                    
      f=dsqrt(dble((lp1-m)*(l+m)))
      con=-con/f
      x(mm)=x(mm)*con
  230 x1(mm)=-m*cot*x(mm)-f*x(m)
      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
      if(l.gt.0) x1(2)=-0.5d0*con*dsqrt(dble(l*lp1))
      return
      end      

      subroutine readmod(io)
c  read in aspherical model
c  note that there must be 40 layers or less
c  and lmax must be .le. 63  and nfun .le. 40
c  the maximum spline order is 40
c  dimension of cof is cof(2*nfun,(lmax+1)(lmax+2)/2,nlay)
      implicit real*8(a-h,o-z)
      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)
      character*400 filnam,linestr
      print *,'enter aspherical model filename: '
      read(*,'(a)') filnam
      open(io,file=filnam)
      read(io,*) nlay
      do 10 ilay=1,nlay
        read(io,'(a)') linestr
        call pstr(linestr,radii(1,ilay),iopt(ilay),lmax(ilay),
     +     ipow(ilay),nfun(ilay))
        istop = (lmax(ilay)+1)*(lmax(ilay)+2)/2
        kread=2*nfun(ilay)
        do 20 ind=1,istop
   20     read(io,*) l,m,(cof(k,ind,ilay),k=1,kread)
   10  continue
      close(io)
      return
      end


      subroutine pstr(stin,rvec,i1,i2,i3,i4)

      implicit real*8(a-h,o-z)
      dimension rvec(*)
      character*400 stin

      imax=0
      k=1
 10   i=k
 15   if (i.gt.400) goto 30
	if (stin(i:i).ne.' ') then
	   imax=imax+1
	   j=i+1
 20	   if (j.gt.400) goto 30
	   if (stin(j:j).eq.' ') then
		k=j+1
		goto 10
	   else
		j=j+1
		goto 20
	   endif
        endif
        i=i+1
	goto 15
      
 30   read(stin,*)(rvec(j),j=1,imax-4),i1,i2,i3,i4

      return
      end


      subroutine radfunx(func,r,ilay)
c Returns values of radial functions at specified radius.
c input:
c  iopt = choice of radial functions
c         1 = Legendre polynomials between knots
c         2 = constant layer
c         3 = lin interp between knots
c         4 = Chebyshev polynomials
c	  5 = natural cubic B-splines
c   r    = radius (km)
c  
c output:
c   func() = radial function values
      implicit real*8(a-h,o-z)
      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)
      dimension func(41),con(5),conc(41)
      data con/
     +  0.7071067811865476d0,1.2247448713915889d0,
     +  0.7905694150420949d0,0.9354143466934853d0,
     +  0.2651650429449553d0/
      data conc /0.7071067811865475d0,1.2247448713915889d0,
     +  1.0350983390135313d0,1.0145993123917847d0,
     +  1.0080322575483707d0,1.0050890913907349d0,
     +  1.0035149493261806d0,1.0025740068320432d0,
     +  1.0019665702377578d0,1.0015515913132542d0,
     +  1.0012554932753530d0,1.0010368069140518d0,
     +  1.0008707010791882d0,1.0007415648034326d0,
     +  1.0006391819124996d0,1.0005566379501474d0,
     +  1.0004891171728023d0,1.0004331817400482d0,
     +  1.0003863241403532d0,1.0003466805443029d0,
     +  1.0003128421787779d0,1.0002837281941048d0,
     +  1.0002584981302063d0,1.0002364904845644d0,
     +  1.0002171788493410d0,1.0002001401000729d0,
     +  1.0001850309942748d0,1.0001715707312959d0,
     +  1.0001595277987336d0,1.0001487099420818d0,
     +  1.0001389564378276d0,1.0001301320846161d0,
     +  1.0001221224893115d0,1.0001148303385539d0,
     +  1.0001081724271648d0,1.0001020772727611d0,
     +  1.0000964831880290d0,1.0000913367129822d0,
     +  1.0000865913323700d0,1.0000822064204604d0,1.0000781463682669d0/

  10  do i=1,41
      func(i)=0.d0
      end do


c  fill up func for function type

      if (iopt(ilay).eq.1) then
       if (r.lt.radii(1,ilay).or.r.gt.radii(2,ilay)) return
        u=(r+r-radii(2,ilay)-radii(1,ilay))/
     +              (radii(2,ilay)-radii(1,ilay))
        usq=u*u
        func(1) = con(1)
        func(2) = con(2)*u
        func(3) = con(3)*(3.d0*usq-1.d0)
        func(4) = con(4)*u*(5.d0*usq-3.d0)
        func(5) = con(5)*(usq*(35.d0*usq-30.d0)+3.d0)
        return
      elseif (iopt(ilay).eq.2) then     
       if (r.lt.radii(1,ilay).or.r.gt.radii(2,ilay)) return
        func(1) = 1.d0
        return
      elseif (iopt(ilay).eq.3) then                
       if (r.lt.radii(1,ilay).or.r.gt.radii(2,ilay)) return
        u = (r - radii(1,ilay)) / (radii(2,ilay) - radii(1,ilay))
        func(1) = u
        func(2) = 1.d0 - u
        return
      elseif (iopt(ilay).eq.4) then                
       if (r.lt.radii(1,ilay).or.r.gt.radii(2,ilay)) return
        u=(r+r-radii(2,ilay)-radii(1,ilay))/
     +              (radii(2,ilay)-radii(1,ilay))
        func(1)= 1.d0 
        func(2)= u  
        u2 = 2.d0 * u
        do 30 i=3,41
   30     func(i) =  u2 * func(i-1) - func(i-2)
        do 35 i=1,41
   35     func(i)=conc(i)*func(i)
      elseif (iopt(ilay).eq.5) then
c add a line to clip at the edges of the layer
        i2=nfun(ilay)+2
        if (r.lt.radii(3,ilay).or.r.gt.radii(i2,ilay)) return
	do j=1,nfun(ilay)
	   func(j)=bsplinen(j-3,r,ilay)
	end do
      endif
      return
      end



       real*8 function bsplinen(i,x,ilay)

c  evaluates natural cubic bspline at ordinate x
c  i is the spline translation index, which
c  should be between -2 and N-2 for the natural 
c  spline basis.


        implicit real*8(a-h,o-z)
      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)


c handle end splines

        nf=nfun(ilay)
	if (i.eq.-2) then
          q0=(radii(4,ilay)-radii(1,ilay))/
     +		(radii(5,ilay)-radii(2,ilay))
          bsplinen=bspline(i,x,ilay)+(1.+q0)*bspline(-3,x,ilay)
	  return
	elseif (i.eq.-1) then
          q0=(radii(4,ilay)-radii(1,ilay))/
     +		(radii(5,ilay)-radii(2,ilay))
          bsplinen=bspline(i,x,ilay)-q0*bspline(-3,x,ilay)
	  return
	elseif (i.eq.nf-4) then
          qN=(radii(nf+4,ilay)-radii(nf+1,ilay))/
     +		(radii(nf+3,ilay)-radii(nf,ilay))
          bsplinen=bspline(nf-4,x,ilay)-qN*bspline(nf-2,x,ilay)
	  return
	elseif (i.eq.nf-3) then
          qN=(radii(nf+4,ilay)-radii(nf+1,ilay))/
     +		(radii(nf+3,ilay)-radii(nf,ilay))
          bsplinen=bspline(nf-3,x,ilay)+(1.+qN)*bspline(nf-2,x,ilay)
	  return
	else
	  bsplinen=bspline(i,x,ilay)
	  return
        endif

	end


       real*8 function bspline(i,x,ilay)

c  evaluates cubic bspline at ordinate x
c  i is the spline translation index, which
c  should be between -3 and N-1 for the 
c  B-spline basis.

        implicit real*8(a-h,o-z)

      parameter (maxlay=30,maxl=16,iradmax=30)
      common/vmap/cof(2*iradmax,(maxl+1)*(maxl+2)/2,maxlay),
     +	radii(iradmax+4,maxlay),nlay,lmax(maxlay),
     + ipow(maxlay),nfun(maxlay),iopt(maxlay)

c find proper subinterval containing x

	do k=0,nfun(ilay)-1
           tval=(x-radii(k+4,ilay))/(radii(k+3,ilay)-radii(k+4,ilay))
	   if (tval.le.1.0.and.tval.ge.0.0) goto 10
        end do

c oops, we are off the end of the interval of definition;
c   return a value of zero

        bspline=0.d0
	return

c k is the interval index, now compute the subinterval for
c  the desired spline translation of i; the subinterval
c  ranges from 0 to 3

 10     j=k-i

	if (j.lt.0) then
	   bspline=0.d0
	   return
	elseif (j.eq.0) then
	   rnum1=x-radii(i+3,ilay)
	   den1=radii(i+6,ilay)-radii(i+3,ilay)
	   den2=radii(i+5,ilay)-radii(i+3,ilay)
	   den3=radii(i+4,ilay)-radii(i+3,ilay)
	   bspline=rnum1*rnum1*rnum1/den1/den2/den3
	   return
	elseif (j.eq.1) then
	   rnum1=x-radii(i+3,ilay)
	   rnum2=radii(i+5,ilay)-x
	   rnum3=radii(i+6,ilay)-x
	   rnum4=x-radii(i+4,ilay)
	   rnum5=radii(i+7,ilay)-x
	   den1=radii(i+6,ilay)-radii(i+3,ilay)
	   den2=radii(i+5,ilay)-radii(i+3,ilay)
	   den3=radii(i+5,ilay)-radii(i+4,ilay)
	   den4=radii(i+6,ilay)-radii(i+4,ilay)
	   den5=radii(i+7,ilay)-radii(i+4,ilay)
	   bspline=rnum1/den1*(rnum1*rnum2/den2/den3+
     +	     rnum3*rnum4/den4/den3)+rnum5*rnum4*rnum4
     +	     /den5/den4/den3
	   return
	elseif (j.eq.2) then
	   rnum1=x-radii(i+3,ilay)
	   rnum2=radii(i+5,ilay)-x
	   rnum3=radii(i+6,ilay)-x
	   rnum4=x-radii(i+4,ilay)
	   rnum5=radii(i+7,ilay)-x
	   den1=radii(i+6,ilay)-radii(i+3,ilay)
	   den2=radii(i+5,ilay)-radii(i+3,ilay)
	   den3=radii(i+5,ilay)-radii(i+4,ilay)
	   den4=radii(i+6,ilay)-radii(i+4,ilay)
	   den5=radii(i+7,ilay)-radii(i+4,ilay)
	   den6=radii(i+6,ilay)-radii(i+5,ilay)
	   den7=radii(i+7,ilay)-radii(i+5,ilay)
	   bspline=rnum5/den5*(rnum4*rnum3/den4/den6-
     +	     rnum5*rnum2/den7/den6)+rnum1*rnum3*rnum3
     +	     /den1/den4/den6
	   return
	elseif (j.eq.3) then
	   rnum5=radii(i+7,ilay)-x
	   den5=radii(i+7,ilay)-radii(i+4,ilay)
	   den7=radii(i+7,ilay)-radii(i+5,ilay)
	   den8=radii(i+7,ilay)-radii(i+6,ilay)
	   bspline=rnum5*rnum5*rnum5/den7/den5/den8
	   return	
	else
	   bspline=0.d0
	   return
	endif

	end
	      

