implicit real*8(a-h,o-z)
c  test program
      call readmod(7)
      print *,'enter depth'
      read(*,*) rad
      rad=6371-rad
      call evmap(rad)
      stop
      end

      subroutine evmap(radius)
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)
      common/vmap/cof(82,231,40),radii(40,40),nlay,lmax(40),ipow(40),
     +   nfun(40),iopt(40)
      data pi/3.14159265358979d0/
      open(37,file='outlm')
      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 continue
      call radfunx(func,radius,ilay)
      kopt=ipow(ilay)
      ind=1
      ipar=iopt(ilay)
      do 40 lp1=1,lmax(ilay)+1
        l=lp1-1
        ind=ind+l
        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
        m=0
        ap=0
        write(37,900)l,m,rp,ap
  900 format(2i5,2g14.6)
        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
        write(37,900)l,m,rp,ap
   45   continue
   40 continue
      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)
      common/vmap/cof(82,231,40),radii(40,40),nlay,lmax(40),ipow(40),
     +   nfun(40),iopt(40)
      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   write(6,*)imax

      read(stin,*)(rvec(j),j=1,imax-4),i1,i2,i3,i4

      write(6,*)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)
      common/vmap/cof(82,231,40),radii(40,40),nlay,lmax(40),ipow(40),
     +   nfun(40),iopt(40)
      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)
      common/vmap/cof(82,231,40),radii(40,40),nlay,lmax(40),ipow(40),
     +   nfun(40),iopt(40)

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)
      common/vmap/cof(82,231,40),radii(40,40),nlay,lmax(40),ipow(40),
     +   nfun(40),iopt(40)

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