ccc pskap.f   Time-stamp: <2004-06-03 12:06:08 hamana>

      implicit real (a-h,o-z)
      implicit integer (i-n)

      common /cosmo/ om,ov,crv,scrv
      common /dpower/rn,gamfin,sig8,anrm
      
      external dw

      parameter (ismin=10)
      parameter (ismax=60)
      parameter (dels=0.1)
      dimension s(ismin:ismax),dps(ismin:ismax),pskap(ismin:ismax)
      dimension dpsnl(ismin:ismax),pskapnl(ismin:ismax)

      parameter (ia=100)
      dimension asf(0:ia),cd(0:ia),d(0:ia),dod(0:ia)

      parameter (nrkmin=-40,nrkmax=60)
      parameter (nrk=nrkmax-nrkmin+1)
      parameter (delrk=0.1)
      dimension rk(nrk),pnl(nrk),pnl2(nrk)

      pi=4.0*atan(1.0)
      twopi=2.0*pi
      tps=2.0*pi**2

      cv=3000.0

      igrfrag=-1
      
      write(*,*)'### pskap.f'
      write(*,*)'### version-1 [3/6/2004]'
      write(*,*)'###'
      write(*,*)'### Use this at your own risk' 
      write(*,*)'### Bug report to hamanatk [at-mark] cc.nao.ac.jp'
      write(*,*)'###'
      write(*,*)'### This uses BBKS CDM P(k) with '
      write(*,*)'### Sugiyama (1995) Gamma factor and'
      write(*,*)'### Smith et al (2003) nonlinear fitting function'
      write(*,*)'###'
      write(*,*)'### This will output a data file pskap.dat'
      write(*,*)'###'
      write(*,*)'### Hubble parameter is defined by'
      write(*,*)'### H_0=100h[km/sec/Mpc]'
      write(*,*)'#########################################'
      
 11   write(*,*)'Cosmology 1, 2 or 3 ?'
      write(*,*)'1: Lambda-flat (Omega_m=0.3, Lambda=0.7)'
      write(*,*)'2: Eds         (Omega_m=1.0, Lambda=0.0)'
      write(*,*)'3: something else'
      read(*,*)icosmo
      
      if (icosmo.eq.1) then
         om=0.3
         ov=0.7
      elseif (icosmo.eq.2) then
         om=1.0
         ov=0.0
      elseif (icosmo.eq.3) then
         write(*,*)'Omega_m ?'
         read(*,*)om
         write(*,*)'Lambda ?'
         read(*,*)ov
      else
         goto 11
      endif
      crv=om+ov-1.0
      scrv=sqrt(abs(crv))
      kcv=0
      if (crv.gt.0.001) kcv=1
      if (crv.lt.-0.001) kcv=-1
      
 12   write(*,*)'Hubble parameter 1, 2 or 3 ?'
      write(*,*)'1: h=0.7'
      write(*,*)'2: h=0.5'
      write(*,*)'3: something else'
      read(*,*)icosmo
      
      if (icosmo.eq.1) then
         hubp=0.7
      elseif (icosmo.eq.2) then
         hubp=0.5
      elseif (icosmo.eq.3) then
         write(*,*)'h ?'
         read(*,*)hubp
      else
         goto 12
      endif

 13   write(*,*)'Baryon density 1, 2 or 3 ?'
      write(*,*)'1: Omega_b=0.04'
      write(*,*)'2: Omega_b=0.0'
      write(*,*)'3: something else'
      read(*,*)icosmo
      
      if (icosmo.eq.1) then
         ob=0.04
      elseif (icosmo.eq.2) then
         ob=0.0
      elseif (icosmo.eq.3) then
         write(*,*)'Omega_b ?'
         read(*,*)ob
      else
         goto 13
      endif

ccc gamma factor by Sugiyama (1995)
      gamfin=exp(ob+sqrt(2.0*hubp)*ob/om)/(om*hubp*cv)

 14   write(*,*)'sigma_8 1, 2 or 3 ?'
      write(*,*)'1: sigma_8=0.9'
      write(*,*)'2: sigma_8=0.5'
      write(*,*)'3: something else'
      read(*,*)icosmo
      
      if (icosmo.eq.1) then
         sig8=0.9
      elseif (icosmo.eq.2) then
         sig8=0.5
      elseif (icosmo.eq.3) then
         write(*,*)'sigma_8 ?'
         read(*,*)sig8
      else
         goto 14
      endif

ccccc normalization of power spectrum      
      call normPk
      coef=9.0/4.0*om*om

 15   write(*,*)'source redshift 1, 2 or 3 ?'
      write(*,*)'1: zs=1.0'
      write(*,*)'2: zs=2.0'
      write(*,*)'3: something else'
      read(*,*)icosmo
      
      if (icosmo.eq.1) then
         zmax=1.0
      elseif (icosmo.eq.2) then
         zmax=2.0
      elseif (icosmo.eq.3) then
         write(*,*)'zs ?'
         read(*,*)zmax
      else
         goto 15
      endif

ccc source redshift
      amin=1.0/(1.0+zmax)
      dela=(1.0-amin)/ia
      dela2=0.5*dela

ccc  complute a-w a-r relation 
      asf(0)=1.0
      cd(0)=0.0
      d(0)=0.0
      do j=1,ia
         asf(j)=1.0-j*dela
         call qromb(dw,asf(j),1.0,ss)
         cd(j)=ss
         d(j)=fk(ss,kcv)
      enddo

      wi=cd(ia)
      fks=1.0/d(ia)
      do j=0,ia  
         wij=wi-cd(j)
         dod1=fk(wij,kcv)*fks
         dod(j)=dod1*dod1
      enddo

      do i=ismin,ismax
         s(i)=10.0**(i*dels)
         pskap(i)=0.0
         dps(i)=0.0
         pskapnl(i)=0.0
         dpsnl(i)=0.0
      enddo

      do i=1,ia
         z=1.0/asf(i)-1.0
         gr2=gr(asf(i),igrfrag)**2
         weight=dod(i)/sqrt(e(asf(i)))
         weightnl=weight/asf(i)**2
         
         call halofit_v3(om,ov,sig8,ob,hubp,z,
     $        nrkmin,nrkmax,nrk,delrk,rk,pnl)
         yp1=(pnl(2)-pnl(1))/(rk(2)-rk(1))
         ypn=(pnl(nrk)-pnl(nrk-1))/(rk(nrk)-rk(nrk-1))
         call spline(rk,pnl,nrk,yp1,ypn,pnl2)
         pow=(log(pnl(nrk))-log(pnl(nrk-1)))/
     $        (log(rk(nrk))-log(rk(nrk-1)))
         
         do ks=ismin,ismax
            ak=s(ks)/d(i)
!!! linear
            dps2=weight*gr2*plin(ak)
            pskap(ks)=pskap(ks)+
     $           (dps(ks)+dps2)*dela2
            dps(ks)=dps2
!!! non-linear
            x=ak/cv
            y=spfit(rk,pnl,pnl2,nrk,pow,x)*tps/ak**3
            dps2=weightnl*y
            pskapnl(ks)=pskapnl(ks)+
     $           (dpsnl(ks)+dps2)*dela2
            dpsnl(ks)=dps2
         enddo
      enddo

      open(10,file='pskap.dat',status='unknown')
      write(10,'(a41)')'# s; Plin; s^2*Plin/2pi; Pnl; s^2*Pnl/2pi' 
      do ks=ismin,ismax
         ps2d=coef*pskap(ks)
         delps2d=ps2d*s(ks)**2/twopi
         ps2dnl=coef*pskapnl(ks)
         delps2dnl=ps2dnl*s(ks)**2/twopi
         write(10,'(5(1pe13.5))')s(ks),ps2d,delps2d,ps2dnl,delps2dnl
      enddo
      close(10)

      end

!!!
      function spfit(rk,pnl,pnl2,nrk,pow,x)
      implicit real (a-h,o-z)
      implicit integer (i-n)

      dimension rk(nrk),pnl(nrk),pnl2(nrk)
      
      if (x.lt.rk(1)) then
         y=pnl(1)*(x/rk(1))**4
      elseif (x.gt.rk(nrk)) then
         y=pnl(nrk)*(x/rk(nrk))**pow
      else
         call splint(rk,pnl,pnl2,nrk,x,y)
      endif
      spfit=y
      return
      end
      
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine halofit_v3(om_m0,om_v0,sig80,om_b0,hubp,z,
     $     nrkmin,nrkmax,ndat,delrk,rkdim,pnldim)

      implicit none 

      real rn,plin,pq,ph,pnl,rk,p_cdm,aexp,z
      real rknl,rneff,rncur,d1,d2,rad,sig,rknl_int_pow
      real om_m,om_v,om_b,h,p_index,gams,sig8,amp
      real om_m0,om_v0,omega_m,omega_v,grow,grow0,gg
      real rn_pd,rn_cdm,pnl_pd,rklin,f_pd
      real f1a,f2a,f3a,f4a,f1b,f2b,f3b,f4b,frac,f1,f2,f3,f4
      real diff,xlogr1,xlogr2,rmid
      real om_b0,hubp
      real rkmin,rkmax

      integer i,j,k,ndat

      real sig80,delrk
      real rkdim(ndat),plindim(ndat),pnldim(ndat)
      integer nrkmin,nrkmax

      common/cospar/om_m,om_v,h,p_index,gams,sig8,amp
      
      sig8=sig80
      gams=om_m0*hubp/exp(om_b0+sqrt(2.0*hubp)*om_b0/om_m0)

      aexp=1./(1.+z)            ! expansion factor

c calculate matter density, vacuum density at desired redshift

      om_m=omega_m(aexp,om_m0,om_v0) 
      om_v=omega_v(aexp,om_m0,om_v0)

c calculate the amplitude of the power spectrum at desired redshift 
c using linear growth factors (gg given by Caroll, Press & Turner 1992, ARAA, 30, 499)

      grow=gg(om_m,om_v) 
      grow0=gg(om_m0,om_v0)       
      amp=aexp*grow/grow0

c calculate nonlinear wavenumber (rknl), effective spectral index (rneff) and 
c curvature (rncur) of the power spectrum at the desired redshift, using method 
c described in Smith et al (2002).

!      write(*,*) 'computing effective spectral quantities:'

      xlogr1=-2.0
      xlogr2=3.5
 10   rmid=(xlogr2+xlogr1)/2.0
      rmid=10**rmid
      call wint(rmid,sig,d1,d2)
      diff=sig-1.0
      if (abs(diff).le.0.001) then
         rknl=1./rmid
         rneff=-3-d1
         rncur=-d2                  
      elseif (diff.gt.0.001) then
         xlogr1=log10(rmid)
         goto 10
      elseif (diff.lt.-0.001) then
         xlogr2=log10(rmid)
         goto 10
      endif

!      write(*,20) 'rknl [h/Mpc] =',rknl,'rneff=',rneff, 'rncur=',rncur
! 20   format(a14,f12.6,2x,a6,f12.6,2x,a6,f12.6)

c now calculate power spectra for a logarithmic range of wavenumbers (rk)

!      write(*,*) 'press return to compute halofit nonlinear power'
!      read(*,*)

      rkmin=nrkmin*delrk
      rkmax=nrkmax*delrk
      j=0
      do i=nrkmin,nrkmax
         j=j+1
         rk=10.0**(i*delrk)

c linear power spectrum !! Remeber => plin = k^3 * P(k) * constant
c constant = 4*pi*V/(2*pi)^3 

         plin=amp*amp*p_cdm(rk,gams,sig8)

c calculate nonlinear power according to halofit: pnl = pq + ph,
c where pq represents the quasi-linear (halo-halo) power and 
c where ph is represents the self-correlation halo term. 
 
         call halofit(rk,rneff,rncur,rknl,plin,pnl,pq,ph)   ! halo fitting formula 
         
         rkdim(j)=rk
         plindim(j)=plin
         pnldim(j)=pnl

      enddo
      
      return
      end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c halo model nonlinear fitting formula as described in 
c Appendix C of Smith et al. (2002)

      subroutine halofit(rk,rn,rncur,rknl,plin,pnl,pq,ph)

      implicit none

      real gam,a,amod,b,c,xmu,xnu,alpha,beta,f1,f2,f3,f4
      real rk,rn,plin,pnl,pq,ph
      real om_m,om_v,h,p_index,gams,sig8,amp
      real rknl,y,rncur,p_cdm
      real f1a,f2a,f3a,f4a,f1b,f2b,f3b,f4b,frac
      real yy

      common/cospar/om_m,om_v,h,p_index,gams,sig8,amp

      gam=0.86485+0.2989*rn+0.1631*rncur
      a=1.4861+1.83693*rn+1.67618*rn*rn+0.7940*rn*rn*rn+
     &0.1670756*rn*rn*rn*rn-0.620695*rncur
      a=10**a      
      b=10**(0.9463+0.9466*rn+0.3084*rn*rn-0.940*rncur)
      c=10**(-0.2807+0.6669*rn+0.3214*rn*rn-0.0793*rncur)
      xmu=10**(-3.54419+0.19086*rn)
      xnu=10**(0.95897+1.2857*rn)
      alpha=1.38848+0.3701*rn-0.1452*rn*rn
      beta=0.8291+0.9854*rn+0.3400*rn**2

      if(abs(1-om_m).gt.0.01) then ! omega evolution 
         f1a=om_m**(-0.0732)
         f2a=om_m**(-0.1423)
         f3a=om_m**(0.0725)
         f1b=om_m**(-0.0307)
         f2b=om_m**(-0.0585)
         f3b=om_m**(0.0743)       
         frac=om_v/(1.-om_m) 
         f1=frac*f1b + (1-frac)*f1a
         f2=frac*f2b + (1-frac)*f2a
         f3=frac*f3b + (1-frac)*f3a
      else         
         f1=1.0
         f2=1.
         f3=1.
      endif

      y=(rk/rknl)

      ph=a*y**(f1*3)/(1+b*y**(f2)+(f3*c*y)**(3-gam))
      ph=ph/(1+xmu*y**(-1)+xnu*y**(-2))
      yy=-y/4.0-y**2/8.0
      if (yy.gt.-88.0) then
         pq=plin*(1+plin)**beta/(1+plin*alpha)*exp(-y/4.0-y**2/8.0)
      else
         pq=0.0
      endif
      pnl=pq+ph

      return
      end       

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c Bond & Efstathiou (1984) approximation to the linear CDM power spectrum 
!!! changed to BBKS96 in v2

      function p_cdm(rk,gams,sig8)
      implicit none 
      real p_cdm,rk,gams,sig8,p_index,rkeff,q,q8,tk,tk8
      integer flag
      save flag,q8,tk8
      if (flag.eq.0) then
         rkeff=0.172+0.011*log(gams/0.36)*log(gams/0.36)
         q8=1.e-20 + rkeff/gams
         tk8=log(1.0+2.34*q8)/((2.34*q8*(1.0+3.89*q8+(16.1*q8)**2+
     $        (5.46*q8)**3+(6.71*q8)**4)**0.25))
         flag=1
      endif
      q=1.e-20 + rk/gams
      tk=log(1.0+2.34*q)/(
     $     (2.34*q*
     $     (1.0+3.89*q+(16.1*q)**2+(5.46*q)**3+(6.71*q)**4)**0.25))
      p_cdm=sig8*sig8*((q/q8)**4)*tk*tk/(tk8*tk8)
      return
      end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c The subroutine wint, finds the effective spectral quantities
c rknl, rneff & rncur. This it does by calculating the radius of 
c the Gaussian filter at which the variance is unity = rknl.
c rneff is defined as the first derivative of the variance, calculated 
c at the nonlinear wavenumber and similarly the rncur is the second
c derivative at the nonlinear wavenumber. 

      subroutine wint(r,sig,d1,d2)
      
      implicit none
      real sum1,sum2,sum3,t,y,x,w1,w2,w3
      real rk, p_cdm, r, sig, d1,d2
      real om_m,om_v,h,p_index,gams,sig8,amp
      integer i,nint

      common/cospar/om_m,om_v,h,p_index,gams,sig8,amp

      nint=3000
      sum1=0.
      sum2=0.
      sum3=0.
      do i=1,nint
         t=(float(i)-0.5)/float(nint)
         y=-1.+1./t
         rk=y
         d2=amp*amp*p_cdm(rk,gams,sig8)
         x=y*r
         w1=exp(-x*x)
         w2=2*x*x*w1
         w3=4*x*x*(1-x*x)*w1
         sum1=sum1+w1*d2/y/t/t
         sum2=sum2+w2*d2/y/t/t
         sum3=sum3+w3*d2/y/t/t
      enddo
      sum1=sum1/float(nint)
      sum2=sum2/float(nint)
      sum3=sum3/float(nint)
      sig=sqrt(sum1)
      d1=-sum2/sum1
      d2=-sum2*sum2/sum1/sum1 - sum3/sum1
      
      return
      end
      
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c evolution of omega matter with expansion factor

      function omega_m(aa,om_m0,om_v0)
      implicit none
      real omega_m,omega_t,om_m0,om_v0,aa
      omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*aa*aa+om_m0/aa)
      omega_m=omega_t*om_m0/(om_m0+om_v0*aa*aa*aa)
      return
      end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c evolution of omega lambda with expansion factor

      function omega_v(aa,om_m0,om_v0)      
      implicit none
      real aa,omega_v,om_m0,om_v0,omega_t
      omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*aa*aa+om_m0/aa)
      omega_v=omega_t*om_v0/(om_v0+om_m0/aa/aa/aa)
      return
      end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c growth factor for linear fluctuations 

      function gg(om_m,om_v)        
      implicit none
      real gg,om_m,om_v
      gg=2.5*om_m/(om_m**(4./7.)-om_v+(1.+om_m/2.)*(1.+om_v/70.))
      return
      end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ccc CDM power spectrum
      function plin(w)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /dpower/rn,gamfin,sig8,anrm
      q=w*gamfin
      plin=anrm*trnf2(q)*w
      return
      end

ccc CDM transfer function by BBKS 1986
      function trnf2(q)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      qq=q*q
      trnf2=0.182628*log(1.0+2.34*q)**2/qq
     $     /sqrt(1.0+3.89*q+259.21*qq
     $     +162.77*qq*q+2027.17*qq*qq)
      return
      end

ccc linear growth factor
      function gr(a,igrfrag)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /cosmo/ om,ov,crv,scrv
      save grnorm
      if (igrfrag.eq.-1) then
         et=1.0/e(1.0)
         ot=om*et
         vt=ov*et
         grnorm=0.4*(ot**(0.571428571)-vt+
     $        (1.0+0.5*ot)*(1.0+vt*0.014285714))/ot
         igrfrag=1
      endif
      et=1.0/e(a)
      ot=om*a*et
      vt=ov*et*a**4
      graw=2.5*ot/(ot**0.571428571-vt+
     $     (1.0+0.5*ot)*(1.0+vt*0.014285714))
      gr=graw*grnorm
      return
      end

ccc linear grwoth rate
      function grpd(a)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /cosmo/ om,ov,crv,scrv
      et=1.0/e(a)
      ot=om*a*et
      vt=ov*et*a**4
      graw=2.5*ot/(ot**0.571428571-vt+
     $     (1.0+0.5*ot)*(1.0+vt*0.014285714))
      grpd=graw
      return
      end

ccc
      function dw(a)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      dw=1.0/sqrt(e(a))
      return
      end

ccc
      function e(a)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /cosmo/ om,ov,crv,scrv
      a2=a*a
      e=a*om-a2*crv+a2*a2*ov
      return
      end

ccc      
      function fk(w,kcv)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /cosmo/ om,ov,crv,scrv
      if (kcv.eq.-1) then
         fk=sinh(scrv*w)/scrv
      elseif (kcv.eq.0) then
         fk=w
      else 
         fk=sin(scrv*w)/scrv
      endif
      return
      end

ccccccc normalization of power spectrum 
      subroutine normPk
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /dpower/rn,gamfin,sig8,anrm
      external delta

      pi=4.0*atan(1.0)
      coe=4.5/pi/pi
      rn=8.0/3000.0

      wmin=-10.0
      wmax=10.0

      call qromb(delta,wmin,wmax,seki)       
      y=coe*seki
      anrm=sig8*sig8/y
      return
      end

ccccccccc k k^2 t^2 W ccccccccccccccccccccccccccccccccccc      
      function delta(x)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      common /dpower/rn,gamfin,sig8,anrm
      w=exp(x)
      q=w*gamfin
      r=w*rn
      filt=(sin(r)-r*cos(r))/r**3
      tranf=trnf2(q)
      delta=w**4*filt*filt*tranf
      return
      end


cccccc integration pakege cccccccccccccccccccccccccccccccc
      SUBROUTINE qromb(func,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL a,b,func,ss,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint,trapzd
      INTEGER j
      REAL dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j)
        if (j.ge.K) then
          call polint(h(j-KM),s(j-KM),K,0.,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      pause 'too many steps in qromb'
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE trapzd(func,a,b,s,n)
      INTEGER n
      REAL a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)pause 'failure in polint'
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.

cccccc spline fitting ccccccccccccccccccccccccccccccccc
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.