Logo Search packages:      
Sourcecode: r-cran-vgam version File versions  Download package

rgam.f

      subroutine nvhb7f(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk, knot,coef,
     &sz,rjcq9o, n9peut,l6xrjt,fpcb2n, sz6ohy, yc1ezl,hts1gp,mk2vyr,
     &thfyl1,la5dcf)
      implicit logical (a-z)
      integer nfiumb4, nk, yc1ezl, hts1gp(3), mk2vyr, thfyl1, la5dcf
      double precision egoxa3, atqh0o, xs(nfiumb4), ys(nfiumb4), ws(
     &nfiumb4), knot(nk+4), coef(nk), sz(nfiumb4), rjcq9o(nfiumb4), 
     &n9peut, l6xrjt, fpcb2n(3), sz6ohy(1)
      call gzmfi3(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk, knot,coef,sz,
     &rjcq9o, n9peut,hts1gp(1),l6xrjt,hts1gp(2), hts1gp(3), fpcb2n(1),
     &fpcb2n(2),fpcb2n(3), yc1ezl, sz6ohy(1), sz6ohy(nk+1),sz6ohy(2*nk+
     &1),sz6ohy(3*nk+1),sz6ohy(4*nk+1), sz6ohy(5*nk+1),sz6ohy(6*nk+1),
     &sz6ohy(7*nk+1),sz6ohy(8*nk+1), sz6ohy(9*nk+1),sz6ohy(9*nk+mk2vyr*
     &nk+1),sz6ohy(9*nk+2*mk2vyr*nk+1), mk2vyr,thfyl1,la5dcf)
      return
      end
      subroutine gzmfi3(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk, knot,coef,
     &sz,rjcq9o, n9peut,rlhz2a,dwgkz6,ispar, pga6nu, lspar,uspar,fjo2dy,
     & yc1ezl, mheq6i, n7cuql,dvpc8x,hdv8br,cbg5ys, vf1jtn,eh6nly,
     &mvx9at,vbxpg4, rlep7v,lunah2,p2ip, mk2vyr,thfyl1,la5dcf)
      implicit logical (a-z)
      integer nfiumb4,nk, rlhz2a,ispar, yc1ezl, mk2vyr,thfyl1,la5dcf
      integer pga6nu
      double precision egoxa3,atqh0o,xs(nfiumb4),ys(nfiumb4),ws(nfiumb4)
     &, knot(nk+4), coef(nk),sz(nfiumb4),rjcq9o(nfiumb4), n9peut,dwgkz6,
     &lspar,uspar,fjo2dy, mheq6i(nk), n7cuql(nk),dvpc8x(nk),hdv8br(nk),
     &cbg5ys(nk), vf1jtn(nk),eh6nly(nk),mvx9at(nk),vbxpg4(nk), rlep7v(
     &mk2vyr,nk),lunah2(mk2vyr,nk),p2ip(thfyl1,nk)
      double precision t1,t2,dyb3po, a,b,c,d,e,kqoy6w,xm,p,q,r,fjo2dy1,
     &fjo2dy2,u,v,w, fu,fv,fw,fx,x, ax,bx
      integer w3gohz, vucgi1r
      double precision hz0fmy, epx9jf
      hz0fmy = 8.0d88
      epx9jf = 0.0d0
      d = 0.5d0
      u = 0.5d0
      dyb3po = 0.5d0
      w3gohz = 1
23000 if(.not.(w3gohz.le.nfiumb4))goto 23002
      if(.not.(ws(w3gohz).gt.0.0d0))goto 23003
      ws(w3gohz) = dsqrt(ws(w3gohz))
23003 continue
       w3gohz = w3gohz+1
      goto 23000
23002 continue
      if(.not.(yc1ezl .eq. 0))goto 23005
      call poqy8c(vf1jtn,eh6nly,mvx9at,vbxpg4,knot,nk)
      call ak9vxi(xs,ys,ws,knot, nfiumb4,nk, mheq6i,n7cuql,dvpc8x,
     &hdv8br,cbg5ys)
      t1 = 0.0d0 
      t2 = 0.0d0
      do 23007 w3gohz = 3,nk-3 
      t1 = t1 + n7cuql(w3gohz) 
23007 continue
      do 23009 w3gohz = 3,nk-3 
      t2 = t2 + vf1jtn(w3gohz) 
23009 continue
      dyb3po = t1/t2
      yc1ezl = 1
23005 continue
      if(.not.(ispar .eq. 1))goto 23011
      call oipu6h(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk,rlhz2a, knot,coef,
     &sz,rjcq9o,n9peut, dwgkz6, mheq6i, n7cuql,dvpc8x,hdv8br,cbg5ys, 
     &vf1jtn,eh6nly,mvx9at,vbxpg4, rlep7v,lunah2,p2ip,mk2vyr,thfyl1,
     &la5dcf)
      return
23011 continue
      ax = lspar 
      bx = uspar
      c = 0.381966011250105097d0
      kqoy6w = 2.0d-5
      vucgi1r = 0
      a = ax
      b = bx
      v = a + c*(b - a)
      w = v
      x = v
      e = 0.0d0
      dwgkz6 = dyb3po * dexp((-2.0d0 + x*6.0d0) * dlog(16.0d0))
      call oipu6h(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk,rlhz2a, knot,coef,
     &sz,rjcq9o,n9peut, dwgkz6, mheq6i, n7cuql,dvpc8x,hdv8br,cbg5ys, 
     &vf1jtn,eh6nly,mvx9at,vbxpg4, rlep7v,lunah2,p2ip,mk2vyr,thfyl1,
     &la5dcf)
      fx = n9peut
      fv = fx
      fw = fx
23013 if(.not.(la5dcf .eq. 0))goto 23014
      vucgi1r = vucgi1r + 1
      xm = 0.5d0*(a + b)
      fjo2dy1 = kqoy6w*dabs(x) + fjo2dy/3.0d0
      fjo2dy2 = 2.0d0*fjo2dy1
      if(.not.((dabs(x - xm) .le. (fjo2dy2 - 0.5d0*(b - a))) .or.(
     &vucgi1r .gt. pga6nu)))goto 23015
      go to 90
23015 continue
      if(.not.((dabs(e) .le. fjo2dy1) .or.(fx .ge. hz0fmy) .or.(fv .ge. 
     &hz0fmy) .or.(fw .ge. hz0fmy)))goto 23017
      go to 40
23017 continue
      r = (x - w)*(fx - fv)
      q = (x - v)*(fx - fw)
      p = (x - v)*q - (x - w)*r
      q = 2.0d0 * (q - r)
      if(.not.(q .gt. 0.0d0))goto 23019
      p = -p
23019 continue
      q = dabs(q)
      r = e
      e = d
      if(.not.((dabs(p) .ge. dabs(0.5d0*q*r)) .or.(q .eq. 0.0d0)))goto 2
     &3021
      go to 40
23021 continue
      if(.not.((p .le. q*(a - x)) .or. (p .ge. q*(b - x))))goto 23023
      go to 40
23023 continue
      d = p/q
      u = x + d
      if(.not.((u - a) .lt. fjo2dy2))goto 23025
      d = dsign(fjo2dy1, xm - x)
23025 continue
      if(.not.((b - u) .lt. fjo2dy2))goto 23027
      d = dsign(fjo2dy1, xm - x)
23027 continue
      go to 50
40    if(.not.(x .ge. xm))goto 23029
      e = a - x
      goto 23030
23029 continue
      e = b - x
23030 continue
      d = c*e
50    if(.not.(dabs(d) .ge. fjo2dy1))goto 23031
      u = x + d
      goto 23032
23031 continue
      u = x + dsign(fjo2dy1, d)
23032 continue
      dwgkz6 = dyb3po * dexp((-2.0d0 + u*6.0) * dlog(16.0d0))
      call oipu6h(egoxa3,atqh0o,xs,ys,ws, nfiumb4,nk,rlhz2a, knot,coef,
     &sz,rjcq9o,n9peut, dwgkz6, mheq6i, n7cuql,dvpc8x,hdv8br,cbg5ys, 
     &vf1jtn,eh6nly,mvx9at,vbxpg4, rlep7v,lunah2,p2ip,mk2vyr,thfyl1,
     &la5dcf)
      fu = n9peut
      if(.not.(fu .gt. hz0fmy))goto 23033
      fu = 2.0d0 * hz0fmy
23033 continue
      if(.not.(fu .le. fx))goto 23035
      if(.not.(u .ge. x))goto 23037
      a = x
      goto 23038
23037 continue
      b = x
23038 continue
      v = w
      fv = fw
      w = x
      fw = fx
      x = u
      fx = fu
      goto 23036
23035 continue
      if(.not.(u .lt. x))goto 23039
      a = u
      goto 23040
23039 continue
      b = u
23040 continue
      if(.not.((fu .le. fw) .or. (w .eq. x)))goto 23041
      v = w
      fv = fw
      w = u
      fw = fu
      goto 23042
23041 continue
      if(.not.((fu .le. fv) .or. (v .eq. x) .or. (v .eq. w)))goto 23043
      v = u
      fv = fu
23043 continue
23042 continue
23036 continue
      goto 23013
23014 continue
90    epx9jf = 0.0d0
      dwgkz6 = dyb3po * dexp((-2.0d0 + x*6.0d0) * dlog(16.0d0))
      n9peut = fx
      return
      return
      end
      subroutine poqy8c(vf1jtn,eh6nly,mvx9at,vbxpg4,tb,nb)
      implicit logical (a-z)
      integer nb
      double precision vf1jtn(nb),eh6nly(nb),mvx9at(nb),vbxpg4(nb),tb(
     &nb+4)
      integer m5xudf,ilo,i6ndbu, ynmzp6, def4wn, nbp1
      integer w3gohz,p1rifj,d9rjek
      integer tlpr2hal
      double precision uq9jtc(4,3),bgu6fw(16),avoe4y(4),yw2(4), wpt
      double precision uoqx2m
      uoqx2m = 1.0d0 / 3.0d0
      ynmzp6 = 3
      def4wn = 4
      nbp1 = nb + 1
      do 23045 w3gohz = 1,nb 
      vf1jtn(w3gohz) = 0.0d0
      eh6nly(w3gohz) = 0.0d0
      mvx9at(w3gohz) = 0.0d0
      vbxpg4(w3gohz) = 0.0d0 
23045 continue
      ilo = 1
      do 23047 w3gohz = 1,nb 
      call vinterv(tb(1), nbp1 ,tb(w3gohz),m5xudf,i6ndbu)
      call vbsplvd(tb,def4wn,tb(w3gohz),m5xudf,bgu6fw,uq9jtc,ynmzp6)
      do 23049 p1rifj = 1,4 
      avoe4y(p1rifj) = uq9jtc(p1rifj,3) 
23049 continue
      call vbsplvd(tb,def4wn,tb(w3gohz+1),m5xudf,bgu6fw,uq9jtc,ynmzp6)
      do 23051 p1rifj = 1,4 
      yw2(p1rifj) = uq9jtc(p1rifj,3) - avoe4y(p1rifj) 
23051 continue
      wpt = tb(w3gohz+1) - tb(w3gohz)
      if(.not.(m5xudf .ge. 4))goto 23053
      do 23055 p1rifj = 1,4 
      d9rjek = p1rifj
      tlpr2hal = m5xudf-4+p1rifj
      vf1jtn(tlpr2hal) = vf1jtn(tlpr2hal) + wpt * (avoe4y(p1rifj)*
     &avoe4y(d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(
     &p1rifj))*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
      d9rjek = p1rifj+1
      if(.not.(d9rjek .le. 4))goto 23057
      eh6nly(tlpr2hal) = eh6nly(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23057 continue
      d9rjek = p1rifj+2
      if(.not.(d9rjek .le. 4))goto 23059
      mvx9at(tlpr2hal) = mvx9at(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23059 continue
      d9rjek = p1rifj+3
      if(.not.(d9rjek .le. 4))goto 23061
      vbxpg4(tlpr2hal) = vbxpg4(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23061 continue
23055 continue
      goto 23054
23053 continue
      if(.not.(m5xudf .eq. 3))goto 23063
      do 23065 p1rifj = 1,3 
      d9rjek = p1rifj
      tlpr2hal = m5xudf-3+p1rifj
      vf1jtn(tlpr2hal) = vf1jtn(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
      d9rjek = p1rifj+1
      if(.not.(d9rjek .le. 3))goto 23067
      eh6nly(tlpr2hal) = eh6nly(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23067 continue
      d9rjek = p1rifj+2
      if(.not.(d9rjek .le. 3))goto 23069
      mvx9at(tlpr2hal) = mvx9at(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23069 continue
23065 continue
      goto 23064
23063 continue
      if(.not.(m5xudf .eq. 2))goto 23071
      do 23073 p1rifj = 1,2 
      d9rjek = p1rifj
      tlpr2hal = m5xudf-2+p1rifj
      vf1jtn(tlpr2hal) = vf1jtn(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
      d9rjek = p1rifj+1
      if(.not.(d9rjek .le. 2))goto 23075
      eh6nly(tlpr2hal) = eh6nly(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23075 continue
23073 continue
      goto 23072
23071 continue
      if(.not.(m5xudf .eq. 1))goto 23077
      do 23079 p1rifj = 1,1 
      d9rjek = p1rifj
      tlpr2hal = m5xudf-1+p1rifj
      vf1jtn(tlpr2hal) = vf1jtn(tlpr2hal) + wpt* (avoe4y(p1rifj)*avoe4y(
     &d9rjek) + (yw2(p1rifj)*avoe4y(d9rjek) + yw2(d9rjek)*avoe4y(p1rifj)
     &)*0.50 + yw2(p1rifj)*yw2(d9rjek)*uoqx2m)
23079 continue
23077 continue
23072 continue
23064 continue
23054 continue
23047 continue
      return
      end
      subroutine gayot2(rlep7v,lunah2,p2ip, mk2vyr,nk,thfyl1,isbkvx6)
      implicit logical (a-z)
      integer mk2vyr,nk,thfyl1,isbkvx6
      double precision rlep7v(mk2vyr,nk), lunah2(mk2vyr,nk), p2ip(
     &thfyl1,nk)
      integer w3gohz, d9rjek, nd6mep
      double precision yrbij3(3),vef2gk(2),cfko0l(1),c0,c1,c2,c3
      double precision wxj6p6, k6nvd6, s6w6ny, ijk1l1, ya6c6v, vj6e6b, 
     &rm44is, pe0ko0, by99io
      c1 = 0.0d0
      c2 = 0.0d0
      c3 = 0.0d0
      yrbij3(1) = 0.0d0
      yrbij3(2) = 0.0d0
      yrbij3(3) = 0.0d0
      vef2gk(1) = 0.0d0
      vef2gk(2) = 0.0d0
      cfko0l(1) = 0.0d0
      do 23081 w3gohz = 1,nk 
      d9rjek = nk-w3gohz+1
      c0 = 1.0d0 / rlep7v(4,d9rjek)
      if(.not.(d9rjek .le. (nk-3)))goto 23083
      c1 = rlep7v(1,d9rjek+3)*c0
      c2 = rlep7v(2,d9rjek+2)*c0
      c3 = rlep7v(3,d9rjek+1)*c0
      goto 23084
23083 continue
      if(.not.(d9rjek .eq. (nk-2)))goto 23085
      c1 = 0.0d0
      c2 = rlep7v(2,d9rjek+2)*c0
      c3 = rlep7v(3,d9rjek+1)*c0
      goto 23086
23085 continue
      if(.not.(d9rjek .eq. (nk-1)))goto 23087
      c1 = 0.0d0
      c2 = 0.0d0
      c3 = rlep7v(3,d9rjek+1)*c0
      goto 23088
23087 continue
      if(.not.(d9rjek .eq. nk))goto 23089
      c1 = 0.0d0
      c2 = 0.0d0
      c3 = 0.0d0
23089 continue
23088 continue
23086 continue
23084 continue
      wxj6p6 = c1*yrbij3(1)
      k6nvd6 = c2*yrbij3(2)
      s6w6ny = c3*yrbij3(3)
      ijk1l1 = c1*yrbij3(2)
      ya6c6v = c2*vef2gk(1)
      vj6e6b = c3*vef2gk(2)
      rm44is = c1*yrbij3(3)
      pe0ko0 = c2*vef2gk(2)
      by99io = c3*cfko0l(1)
      lunah2(1,d9rjek) = 0.0d0 - (wxj6p6+k6nvd6+s6w6ny)
      lunah2(2,d9rjek) = 0.0d0 - (ijk1l1+ya6c6v+vj6e6b)
      lunah2(3,d9rjek) = 0.0d0 - (rm44is+pe0ko0+by99io)
      lunah2(4,d9rjek) = c0**2 + c1*(wxj6p6 + 2.0d0*(k6nvd6 + s6w6ny)) +
     & c2*(ya6c6v + 2.0d0* vj6e6b) + c3*by99io
      yrbij3(1) = vef2gk(1)
      yrbij3(2) = vef2gk(2)
      yrbij3(3) = lunah2(2,d9rjek)
      vef2gk(1) = cfko0l(1)
      vef2gk(2) = lunah2(3,d9rjek)
      cfko0l(1) = lunah2(4,d9rjek)
23081 continue
      if(.not.(isbkvx6 .eq. 0))goto 23091
      return
23091 continue
      do 23093 w3gohz = 1,nk 
      d9rjek = nk-w3gohz+1
      nd6mep = 1
23095 if(.not.(nd6mep.le.4.and.d9rjek+nd6mep-1.le.nk))goto 23097
      p2ip(d9rjek,d9rjek+nd6mep-1) = lunah2(5-nd6mep,d9rjek)
       nd6mep = nd6mep+1
      goto 23095
23097 continue
23093 continue
      do 23098 w3gohz = 1,nk 
      d9rjek = nk-w3gohz+1
      nd6mep = d9rjek-4
23100 if(.not.(nd6mep.ge.1))goto 23102
      c0 = 1.0 / rlep7v(4,nd6mep) 
      c1 = rlep7v(1,nd6mep+3)*c0
      c2 = rlep7v(2,nd6mep+2)*c0 
      c3 = rlep7v(3,nd6mep+1)*c0
      p2ip(nd6mep,d9rjek) = 0.0d0- ( c1*p2ip(nd6mep+3,d9rjek) + c2*p2ip(
     &nd6mep+2,d9rjek) + c3*p2ip(nd6mep+1,d9rjek) )
       nd6mep = nd6mep-1
      goto 23100
23102 continue
23098 continue
      return
      end
      subroutine oipu6h(egoxa3,atqh0o,x,y,w, nfiumb4,nk,rlhz2a, knot,
     &coef,sz,rjcq9o, n9peut, dwgkz6, mheq6i, n7cuql,dvpc8x,hdv8br,
     &cbg5ys, vf1jtn,eh6nly,mvx9at,vbxpg4, rlep7v,lunah2,p2ip, mk2vyr,
     &thfyl1,fjg0qv)
      implicit logical (a-z)
      integer nfiumb4,nk,rlhz2a, mk2vyr,thfyl1,fjg0qv
      double precision egoxa3,atqh0o,x(nfiumb4),y(nfiumb4),w(nfiumb4)
      double precision knot(nk+4), coef(nk),sz(nfiumb4),rjcq9o(nfiumb4),
     & n9peut, dwgkz6, mheq6i(nk)
      double precision n7cuql(nk),dvpc8x(nk),hdv8br(nk),cbg5ys(nk)
      double precision vf1jtn(nk),eh6nly(nk),mvx9at(nk),vbxpg4(nk), 
     &rlep7v(mk2vyr,nk),lunah2(mk2vyr,nk),p2ip(thfyl1,nk)
      double precision das4bx, bgu6fw(16), b0,b1,b2,b3,kqoy6w, uq9jtc(4,
     &1), xv,eqdf
      double precision zh0bs0
      double precision risyv0
      integer oht3ga, ynmzp6, ilo, i6ndbu, d9rjek, w3gohz
      integer px1yhr, m5xudf, def4wn, uxzze7, nkp1
      ilo = 1
      kqoy6w = 0.1d-10
      oht3ga = 0
      ynmzp6 = 3
      def4wn = 4
      uxzze7 = 1
      nkp1 = nk + 1
      do 23103 w3gohz = 1,nk 
      coef(w3gohz) = mheq6i(w3gohz) 
23103 continue
      do 23105 w3gohz = 1,nk 
      rlep7v(4,w3gohz) = n7cuql(w3gohz)+dwgkz6*vf1jtn(w3gohz) 
23105 continue
      do 23107 w3gohz = 1,(nk-1) 
      rlep7v(3,w3gohz+1) = dvpc8x(w3gohz)+dwgkz6*eh6nly(w3gohz) 
23107 continue
      do 23109 w3gohz = 1,(nk-2) 
      rlep7v(2,w3gohz+2) = hdv8br(w3gohz)+dwgkz6*mvx9at(w3gohz) 
23109 continue
      do 23111 w3gohz = 1,(nk-3) 
      rlep7v(1,w3gohz+3) = cbg5ys(w3gohz)+dwgkz6*vbxpg4(w3gohz) 
23111 continue
      call dpbfa8(rlep7v,mk2vyr,nk,ynmzp6,fjg0qv)
      if(.not.(fjg0qv .ne. 0))goto 23113
      return
23113 continue
      call dpbsl8(rlep7v,mk2vyr,nk,ynmzp6,coef)
      px1yhr = 1
      do 23115 w3gohz = 1,nfiumb4 
      xv = x(w3gohz)
      call wbvalue(knot,coef, nk,def4wn,xv,oht3ga, sz(w3gohz))
23115 continue
      if(.not.(rlhz2a .eq. 0))goto 23117
      return
23117 continue
      call gayot2(rlep7v,lunah2,p2ip, mk2vyr,nk,thfyl1,oht3ga)
      do 23119 w3gohz = 1,nfiumb4 
      xv = x(w3gohz)
      call vinterv(knot(1), nkp1 ,xv,m5xudf,i6ndbu)
      if(.not.(i6ndbu .eq. -1))goto 23121
      m5xudf = 4 
      xv = knot(4) + kqoy6w 
23121 continue
      if(.not.(i6ndbu .eq. 1))goto 23123
      m5xudf = nk 
      xv = knot(nk+1) - kqoy6w 
23123 continue
      d9rjek = m5xudf-3
      call vbsplvd(knot,def4wn,xv,m5xudf,bgu6fw,uq9jtc,uxzze7)
      b0 = uq9jtc(1,1)
      b1 = uq9jtc(2,1)
      b2 = uq9jtc(3,1)
      b3 = uq9jtc(4,1)
      zh0bs0 = (b0 *(lunah2(4,d9rjek)*b0 + 2.0d0*(lunah2(3,d9rjek)*b1 + 
     &lunah2(2,d9rjek)*b2 + lunah2(1,d9rjek)*b3)) + b1 *(lunah2(4,
     &d9rjek+1)*b1 + 2.0d0*(lunah2(3,d9rjek+1)*b2 + lunah2(2,d9rjek+1)*
     &b3)) + b2 *(lunah2(4,d9rjek+2)*b2 + 2.0d0* lunah2(3,d9rjek+2)*b3 )
     &+ b3**2* lunah2(4,d9rjek+3)) * w(w3gohz)**2
      rjcq9o(w3gohz) = zh0bs0
23119 continue
      if(.not.(rlhz2a .eq. 1))goto 23125
      das4bx = 0.0d0 
      eqdf = 0.0d0 
      risyv0 = 0.0d0
      do 23127 w3gohz = 1,nfiumb4 
      das4bx = das4bx + ((y(w3gohz)-sz(w3gohz))*w(w3gohz))**2
      eqdf = eqdf + rjcq9o(w3gohz)
      risyv0 = risyv0 + w(w3gohz)*w(w3gohz)
23127 continue
      n9peut = (das4bx/risyv0)/((1.0d0-(atqh0o+egoxa3*eqdf)/risyv0)**2)
      goto 23126
23125 continue
      if(.not.(rlhz2a .eq. 2))goto 23129
      n9peut = 0.0d0
      risyv0 = 0.0d0
      do 23131 w3gohz = 1,nfiumb4 
      n9peut = n9peut + (((y(w3gohz)-sz(w3gohz))*w(w3gohz))/(1.0d0-
     &rjcq9o(w3gohz)))**2
      risyv0 = risyv0 + w(w3gohz)*w(w3gohz)
23131 continue
      n9peut = n9peut / risyv0
      goto 23130
23129 continue
      n9peut = 0.0d0
      do 23133 w3gohz = 1,nfiumb4 
      n9peut = n9peut+rjcq9o(w3gohz)
23133 continue
      n9peut = 3.0d0 + (atqh0o-n9peut)**2
23130 continue
23126 continue
      return
      end
      subroutine ak9vxi(p3vlea,hr83e,w,onyz6j, xl6qgm,nfiumb4, wevr5o,
     &n7cuql,dvpc8x,hdv8br,cbg5ys)
      implicit logical (a-z)
      integer xl6qgm,nfiumb4
      double precision p3vlea(xl6qgm),hr83e(xl6qgm),w(xl6qgm),onyz6j(
     &nfiumb4+4), wevr5o(nfiumb4), n7cuql(nfiumb4),dvpc8x(nfiumb4),
     &hdv8br(nfiumb4),cbg5ys(nfiumb4)
      double precision kqoy6w,uq9jtc(4,1),bgu6fw(16)
      double precision gyn0o0, sce5d5
      integer d9rjek,w3gohz,ilo,m5xudf,i6ndbu, nhwi2tb1
      integer def4wn, uxzze7
      uxzze7 = 1
      def4wn = 4
      nhwi2tb1 = nfiumb4 + 1
      do 23135 w3gohz = 1,nfiumb4 
      wevr5o(w3gohz) = 0.0d0 
      n7cuql(w3gohz) = 0.0d0 
      dvpc8x(w3gohz) = 0.0d0
      hdv8br(w3gohz) = 0.0d0 
      cbg5ys(w3gohz) = 0.0d0
23135 continue
      ilo = 1
      kqoy6w = 0.1d-9
      do 23137 w3gohz = 1,xl6qgm 
      call vinterv(onyz6j(1), nhwi2tb1 ,p3vlea(w3gohz),m5xudf,i6ndbu)
      if(.not.(i6ndbu .eq. 1))goto 23139
      if(.not.(p3vlea(w3gohz) .le. (onyz6j(m5xudf)+kqoy6w)))goto 23141
      m5xudf = m5xudf-1
      goto 23142
23141 continue
      return
23142 continue
23139 continue
      call vbsplvd(onyz6j,def4wn,p3vlea(w3gohz),m5xudf,bgu6fw,uq9jtc,
     &uxzze7)
      d9rjek = m5xudf-4+1
      gyn0o0 = w(w3gohz)**2
      sce5d5 = gyn0o0 * uq9jtc(1,1)
      wevr5o(d9rjek) = wevr5o(d9rjek) + sce5d5*hr83e(w3gohz)
      n7cuql(d9rjek) = n7cuql(d9rjek) + sce5d5*uq9jtc(1,1)
      dvpc8x(d9rjek) = dvpc8x(d9rjek) + sce5d5*uq9jtc(2,1)
      hdv8br(d9rjek) = hdv8br(d9rjek) + sce5d5*uq9jtc(3,1)
      cbg5ys(d9rjek) = cbg5ys(d9rjek) + sce5d5*uq9jtc(4,1)
      d9rjek = m5xudf-4+2
      sce5d5 = gyn0o0 * uq9jtc(2,1)
      wevr5o(d9rjek) = wevr5o(d9rjek) + sce5d5*hr83e(w3gohz)
      n7cuql(d9rjek) = n7cuql(d9rjek) + sce5d5*uq9jtc(2,1)
      dvpc8x(d9rjek) = dvpc8x(d9rjek) + sce5d5*uq9jtc(3,1)
      hdv8br(d9rjek) = hdv8br(d9rjek) + sce5d5*uq9jtc(4,1)
      d9rjek = m5xudf-4+3
      sce5d5 = gyn0o0 * uq9jtc(3,1)
      wevr5o(d9rjek) = wevr5o(d9rjek) + sce5d5*hr83e(w3gohz)
      n7cuql(d9rjek) = n7cuql(d9rjek) + sce5d5*uq9jtc(3,1)
      dvpc8x(d9rjek) = dvpc8x(d9rjek) + sce5d5*uq9jtc(4,1)
      d9rjek = m5xudf
      sce5d5 = gyn0o0 * uq9jtc(4,1)
      wevr5o(d9rjek) = wevr5o(d9rjek) + sce5d5*hr83e(w3gohz)
      n7cuql(d9rjek) = n7cuql(d9rjek) + sce5d5*uq9jtc(4,1)
23137 continue
      return
      end

Generated by  Doxygen 1.6.0   Back to index