45 lines
1.2 KiB
Fortran
45 lines
1.2 KiB
Fortran
subroutine intrp(wltab,absop,wlgrid,abgrd,nfr,nfgrid)
|
|
c =====================================================
|
|
c
|
|
c a more efficient interpolation routine - using bisection
|
|
c
|
|
INCLUDE 'PARAMS.FOR'
|
|
dimension wltab(1),absop(1),wlgrid(1),abgrd(1)
|
|
dimension yint(mfgrid),jint(mfgrid)
|
|
c
|
|
c set up interpolation coefficients for an interpolation
|
|
c by bisection
|
|
c
|
|
fr1=wltab(1)
|
|
fr2=wltab(nfr)
|
|
do ij=1,nfgrid
|
|
xint=wlgrid(ij)
|
|
jl=0
|
|
ju=nfr+1
|
|
10 continue
|
|
if(ju-jl.gt.1) then
|
|
jm=(ju+jl)/2
|
|
if((fr2.gt.fr1).eqv.(xint.gt.wltab(jm))) then
|
|
jl=jm
|
|
else
|
|
ju=jm
|
|
end if
|
|
go to 10
|
|
end if
|
|
j=jl
|
|
if(j.eq.nfr) j=j-1
|
|
if(j.eq.0) j=j+1
|
|
jint(ij)=j
|
|
c yint(ij)=un/log10(wltab(j+1)/wltab(j))
|
|
yint(ij)=1./(wltab(j+1)-wltab(j))
|
|
end do
|
|
c
|
|
do ij=1,nfgrid
|
|
j=jint(ij)
|
|
rc=(absop(j+1)-absop(j))*yint(ij)
|
|
c abgrd(ij)=rc*log10(wlgrid(ij)/wltab(j))+absop(j)
|
|
abgrd(ij)=rc*(wlgrid(ij)-wltab(j))+absop(j)
|
|
end do
|
|
return
|
|
end
|