text:
c
c     Laplace interactions
c
c     real-valued charge, dipstr, pot, grad, hess
c     real-valued dipvec
c
c**********************************************************************
      subroutine r2d_directcp(nd,sources,ns,charge,
     $           targ,nt,pot,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,targ,nt,thresh
cf2py  intent(out) pot
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j log(|targ(*,i)-source(*,j)|)*charge(ii,j)
c
c     The potential is not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)   (real *8)      : potential is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt),xdiff,ydiff,rr,r
      real *8 thresh,rtmp,thresh2
      real *8 pot(nd,nt)
      real *8 charge(nd,ns)

      thresh2 = thresh*thresh
c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            rr=xdiff*xdiff+ydiff*ydiff
            if(rr.le.thresh2) goto 1000

            rtmp = log(rr)/2
            do ii=1,nd
               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
            enddo
 1000       continue
         enddo
      enddo
      return
      end
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directcg(nd,sources,ns,charge,targ,nt,pot,
     1             grad,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,targ,nt,thresh
cf2py  intent(out) pot,grad
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT and gradients GRAD
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c
c     pot(ii,i)  = \sum_j log(|targ(*,i)-source(*,j)|)*charge(ii,j)
c
c     grad(ii,1:2,i)  = \nabla pot(ii,i)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential and gradient
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)      : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh,rtmp,thresh2
      real *8 pot(nd,nt),grad(nd,2,nt)
      real *8 dx,dy
      real *8 charge(nd,ns)

      thresh2 = thresh*thresh
c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            rr=xdiff*xdiff+ydiff*ydiff
            if(rr.lt.thresh2) goto 1000
            rtmp = log(rr)/2
            dx = xdiff/rr
            dy = ydiff/rr
c
            do ii = 1,nd
               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
               grad(ii,1,j) = grad(ii,1,j) + dx*charge(ii,i)
               grad(ii,2,j) = grad(ii,2,j) + dy*charge(ii,i)
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directch(nd,sources,ns,charge,targ,nt,
     1           pot,grad,hess,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,targ,nt,thresh
cf2py  intent(out) pot,grad,hess
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT, gradients GRAD
c     and Hessians at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c
c     pot(ii,i)  = \sum_j log(|targ(*,i)-source(*,j)|)*charge(ii,j)
c
c     grad  = d/dz(pot)
c     hess = d^2/dz^2(pot)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential,gradient,and
c                       hessian
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)      : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c     hess(nd,3,nt)  (real *8)      : Hessian is incremented
c---------------------------------------------------------------------
      integer i,ns,ifexpon,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh,thresh2
      real *8 pot(nd,nt),grad(nd,2,nt),hess(nd,3,nt)
      real *8 rtmp, xdiff2, ydiff2
      real *8 rr2, dx, dy, dxx, dxy, dyy
      real *8 charge(nd,ns)

      thresh2 = thresh*thresh

c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            xdiff2 = xdiff*xdiff
            ydiff2 = ydiff*ydiff
            rr=xdiff2+ydiff2
            rr2 = rr*rr
            if(rr.lt.thresh2) goto 1000
            rtmp = log(rr)/2
            dx = xdiff/rr
            dy = ydiff/rr
            dxx = (rr - 2*xdiff2)/rr2
            dxy = -2*xdiff*ydiff/rr2
            dyy = (rr - 2*ydiff2)/rr2

c
            do ii = 1,nd
               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
               grad(ii,1,j) = grad(ii,1,j) + dx*charge(ii,i)
               grad(ii,2,j) = grad(ii,2,j) + dy*charge(ii,i)
               hess(ii,1,j) = hess(ii,1,j) + dxx*charge(ii,i)
               hess(ii,2,j) = hess(ii,2,j) + dxy*charge(ii,i)
               hess(ii,3,j) = hess(ii,3,j) + dyy*charge(ii,i)
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directdp(nd,sources,ns,dipstr,dipvec,
     $           targ,nt,pot,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                               \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     The potential is not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     dipstr(nd,ns) :   dipole strengths
c     dipvec(nd,2,ns) :   dipole orientations
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)   (real *8)      : potential is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt),xdiff,ydiff,rr,r
      real *8 thresh,thresh2,p1,p2
      real *8 pot(nd,nt)
      real *8 dipstr(nd,ns)
      real *8 dipvec(nd,2,ns)
c
      thresh2 = thresh*thresh

      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            rr = xdiff*xdiff + ydiff*ydiff
            if(rr.le.thresh2) goto 1000

            p1 = -xdiff/rr
            p2 = -ydiff/rr

            do ii=1,nd
               pot(ii,j) = pot(ii,j) + dipstr(ii,i)*(dipvec(ii,1,i)*p1
     1              + dipvec(ii,2,i)*p2)
            enddo
 1000       continue
         enddo
      enddo
      return
      end
c
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directdg(nd,sources,ns,dipstr,dipvec,
     1     targ,nt,pot,grad,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot,grad
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT and gradients GRAD
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                               \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     grad  = \nabla(pot)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     dipstr(nd,ns) :   charge strengths
c     dipvec(nd,2,ns) :   dipole orientations
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential and gradient
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)    : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh,rtmp
      real *8 pot(nd,nt),grad(nd,2,nt)
      real *8 dipstr(nd,ns),d1,d2
      real *8 dx1,dx2,dy1,dy2
      real *8 dipvec(nd,2,ns)
      real *8 xdiff2,ydiff2,p1,p2,thresh2,rr2
c
      thresh2 = thresh*thresh

      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            xdiff2 = xdiff*xdiff
            ydiff2 = ydiff*ydiff
            rr = xdiff2 + ydiff2
            if(rr.le.thresh2) goto 1000

            rr2 = rr*rr

            p1 = -xdiff/rr
            p2 = -ydiff/rr

            dx1 = -(rr-2*xdiff2)/rr2
            dx2 = 2*xdiff*ydiff/rr2
            dy1 = dx2
            dy2 = -(rr-2*ydiff2)/rr2

c
            do ii = 1,nd
               d1 = dipstr(ii,i)*dipvec(ii,1,i)
               d2 = dipstr(ii,i)*dipvec(ii,2,i)
               pot(ii,j) = pot(ii,j) + d1*p1 + d2*p2
               grad(ii,1,j) = grad(ii,1,j) + d1*dx1 + d2*dx2
               grad(ii,2,j) = grad(ii,2,j) + d1*dy1 + d2*dy2
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directdh(nd,sources,ns,dipstr,dipvec,targ,nt,
     1           pot,grad,hess,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot,grad,hess
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT, gradients GRAD
c     and Hessians at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                               \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     grad  = \nabla (pot)
c     hess = \nabla \nabla (pot) (dxx,dxy,dyy)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     dipstr(nd,ns) :   charge strengths
c     dipvec(nd,2,ns) :   dipole orientations
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential,gradient,and
c                       hessian
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)      : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c     hess(nd,3,nt)  (real *8)      : Hessian is incremented
c---------------------------------------------------------------------
      integer i,ns,ifexpon,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh
      real *8 pot(nd,nt),grad(nd,2,nt),hess(nd,3,nt)
      real *8 rtmp
      real *8 dipstr(nd,ns),d1,d2
      real *8 dipvec(nd,2,ns),dx1,dx2,dy1,dy2
      real *8 xdiff2,ydiff2,p1,p2,dxx1,dxx2,dxy1,dxy2,dyy1,dyy2
      real *8 rr2,rr3,thresh2

c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            xdiff2 = xdiff*xdiff
            ydiff2 = ydiff*ydiff
            rr = xdiff2 + ydiff2
            if(rr.le.thresh2) goto 1000

            rr2 = rr*rr
            rr3 = rr2*rr

            p1 = -xdiff/rr
            p2 = -ydiff/rr

            dx1 = -(rr-2*xdiff2)/rr2
            dx2 = 2*xdiff*ydiff/rr2
            dy1 = dx2
            dy2 = -(rr-2*ydiff2)/rr2

            dxx1 = -2*xdiff*(xdiff2-3*ydiff2)/rr3
            dxx2 = 2*ydiff*(ydiff2-3*xdiff2)/rr3
            dxy1 = dxx2
            dxy2 = -dxx1
            dyy1 = dxy2
            dyy2 = -dxx2

c
            do ii = 1,nd
               d1 = dipstr(ii,i)*dipvec(ii,1,i)
               d2 = dipstr(ii,i)*dipvec(ii,2,i)
               pot(ii,j) = pot(ii,j) + d1*p1 + d2*p2
               grad(ii,1,j) = grad(ii,1,j) + d1*dx1 + d2*dx2
               grad(ii,2,j) = grad(ii,2,j) + d1*dy1 + d2*dy2
               hess(ii,1,j) = hess(ii,1,j) + d1*dxx1 + d2*dxx2
               hess(ii,2,j) = hess(ii,2,j) + d1*dxy1 + d2*dxy2
               hess(ii,3,j) = hess(ii,3,j) + d1*dyy1 + d2*dyy2
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directcdp(nd,sources,ns,charge,
     $           dipstr,dipvec,targ,nt,pot,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j log(|z_{j}|)*charge(ii,j)+
c                   dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                      \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     The potential is not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     dipstr(nd,ns) :   dipole strengths
c     dipvec(nd,2,ns) :   dipole orientations
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)   (real *8)      : potential is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt),xdiff,ydiff,rr,r
      real *8 thresh,rtmp,thresh2
      real *8 pot(nd,nt)
      real *8 dipvec(nd,2,ns)
      real *8 charge(nd,ns),dipstr(nd,ns)
      real *8 xdiff2,ydiff2,p1,p2

      thresh2 = thresh*thresh
c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            rr=xdiff*xdiff+ydiff*ydiff
            if(rr.le.thresh2) goto 1000

            rtmp = log(rr)/2.0d0

            p1 = -xdiff/rr
            p2 = -ydiff/rr

            do ii=1,nd
               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
               pot(ii,j) = pot(ii,j) + dipstr(ii,i)*(dipvec(ii,1,i)*p1
     1              + dipvec(ii,2,i)*p2)
            enddo
 1000       continue
         enddo
      enddo
      return
      end
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directcdg(nd,sources,ns,charge,dipstr,dipvec,
     1     targ,nt,pot,grad,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot,grad
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT and gradients GRAD
c     at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c     pot(ii,i)  = \sum_j log(|z_{j}|)*charge(ii,j)+
c                   dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                      \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     grad  = \nabla(pot)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     dipstr(nd,ns) :   dipole strengths
c     dipvec(nd,2,ns) :   dipole orientations
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential and gradient
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)      : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c---------------------------------------------------------------------
      integer i,ns,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh,rtmp,thresh2
      real *8 pot(nd,nt),grad(nd,2,nt)
      real *8 charge(nd,ns),dipstr(nd,ns)
      real *8 d1,d2
      real *8 dipvec(nd,2,ns)
      real *8 xdiff2,ydiff2,p1,p2,dxx1,dxx2,dxy1,dxy2,dyy1,dyy2
      real *8 rr2,rr3,dx,dy,dx1,dx2,dy1,dy2

      thresh2 = thresh*thresh
c
      do j = 1,nt
         do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            xdiff2 = xdiff*xdiff
            ydiff2 = ydiff*ydiff
            rr = xdiff2 + ydiff2
            if(rr.lt.thresh2) goto 1000
            rtmp = log(rr)/2

            rr2 = rr*rr

            dx = xdiff/rr
            dy = ydiff/rr

            p1 = -dx
            p2 = -dy

            dx1 = -(rr-2*xdiff2)/rr2
            dx2 = 2*xdiff*ydiff/rr2
            dy1 = dx2
            dy2 = -(rr-2*ydiff2)/rr2

c
            do ii = 1,nd
               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
               grad(ii,1,j) = grad(ii,1,j) + dx*charge(ii,i)
               grad(ii,2,j) = grad(ii,2,j) + dy*charge(ii,i)

               d1 = dipstr(ii,i)*dipvec(ii,1,i)
               d2 = dipstr(ii,i)*dipvec(ii,2,i)

               pot(ii,j) = pot(ii,j) + d1*p1 + d2*p2
               grad(ii,1,j) = grad(ii,1,j) + d1*dx1 + d2*dx2
               grad(ii,2,j) = grad(ii,2,j) + d1*dy1 + d2*dy2
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c
c
c
c
c
c
c**********************************************************************
      subroutine r2d_directcdh(nd,sources,ns,charge,dipstr,dipvec,
     1     targ,nt,pot,grad,hess,thresh)
cf2py  intent(in) nd
cf2py  intent(in) ns,sources,charge,dipstr,dipvec,targ,nt,thresh
cf2py  intent(out) pot,grad,hess
      implicit none
c**********************************************************************
c
c     This subroutine INCREMENTS the potentials POT, gradients GRAD
c     and Hessians at the target points TARGET, due to a vector of
c     charges at SOURCE(2,ns).
c     We use the unscaled version of log
c     response: i.e., log|z|
c
c
c     pot(ii,i)  = \sum_j log(|z_{j}|)*charge(ii,j)+
c                   dipstr(ii,j)* dipstr(ii,:,j) \cdot
c                      \nabla_src \log |targ(:,i)-sources(:,j)|
c
c     grad  = \nabla(pot)
c     hess = \nabla\nabla (pot)    (dxx,dxy,dyy)
c
c     The potential and gradient are not computed if |r| < thresh
c     (Recommended value for threshold in an FMM is
c     R*eps, where R is the size of the computation
c     domain and eps is machine precision
c
c---------------------------------------------------------------------
c     INPUT:
c
c     sources(2,ns) :   location of the sources
c     ns            :   number of sources
c     charge(nd,ns) :   charge strengths
c     dipstr(nd,ns) :   dipole strengths
c     targ(2,nt)    :   location of the targets
c     thresh        :   threshold for computing potential,gradient,and
c                       hessian
c---------------------------------------------------------------------
c     OUTPUT:
c
c     pot(nd,nt)     (real *8)      : potential is incremented
c     grad(nd,2,nt)  (real *8)      : gradient is incremented
c     hess(nd,3,nt)  (real *8)      : Hessian is incremented
c---------------------------------------------------------------------
      integer i,ns,ifexpon,ii,nd,j,nt
      real *8 sources(2,ns),targ(2,nt)
      real *8 xdiff,ydiff,rr,r,thresh,thresh2
      real *8 pot(nd,nt),grad(nd,2,nt),hess(nd,3,nt)
      real *8 rtmp
      real *8 charge(nd,ns),dipstr(nd,ns),d1,d2
      real *8 dipvec(nd,2,ns),dx1,dx2,dy1,dy2
      real *8 xdiff2,ydiff2,p1,p2,dxx1,dxx2,dxy1,dxy2,dyy1,dyy2
      real *8 rr2,rr3,dxx,dxy,dyy,dx,dy

      thresh2 = thresh*thresh

c
      do j = 1,nt
      do i = 1,ns
            xdiff=targ(1,j)-sources(1,i)
            ydiff=targ(2,j)-sources(2,i)
            xdiff2 = xdiff*xdiff
            ydiff2 = ydiff*ydiff
            rr = xdiff2 + ydiff2
            if(rr.lt.thresh2) goto 1000
            rtmp = log(rr)/2

            rr2 = rr*rr
            rr3 = rr2*rr

            dx = xdiff/rr
            dy = ydiff/rr
            dxx = (rr - 2*xdiff2)/rr2
            dxy = -2*xdiff*ydiff/rr2
            dyy = (rr - 2*ydiff2)/rr2

            p1 = -dx
            p2 = -dy

            dx1 = -dxx
            dx2 = -dxy
            dy1 = dx2
            dy2 = -dyy

            dxx1 = -2*xdiff*(xdiff2-3*ydiff2)/rr3
            dxx2 = 2*ydiff*(ydiff2-3*xdiff2)/rr3
            dxy1 = dxx2
            dxy2 = -dxx1
            dyy1 = dxy2
            dyy2 = -dxx2

            do ii = 1,nd

               pot(ii,j) = pot(ii,j) + rtmp*charge(ii,i)
               grad(ii,1,j) = grad(ii,1,j) + dx*charge(ii,i)
               grad(ii,2,j) = grad(ii,2,j) + dy*charge(ii,i)
               hess(ii,1,j) = hess(ii,1,j) + dxx*charge(ii,i)
               hess(ii,2,j) = hess(ii,2,j) + dxy*charge(ii,i)
               hess(ii,3,j) = hess(ii,3,j) + dyy*charge(ii,i)

               d1 = dipstr(ii,i)*dipvec(ii,1,i)
               d2 = dipstr(ii,i)*dipvec(ii,2,i)

               pot(ii,j) = pot(ii,j) + d1*p1 + d2*p2
               grad(ii,1,j) = grad(ii,1,j) + d1*dx1 + d2*dx2
               grad(ii,2,j) = grad(ii,2,j) + d1*dy1 + d2*dy2
               hess(ii,1,j) = hess(ii,1,j) + d1*dxx1 + d2*dxx2
               hess(ii,2,j) = hess(ii,2,j) + d1*dxy1 + d2*dxy2
               hess(ii,3,j) = hess(ii,3,j) + d1*dyy1 + d2*dyy2
            enddo
 1000       continue
         enddo
      enddo

      return
      end
c
c
c


text: #redirect 3rd & U Cafe

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! ej9a.f90
!
! Copyright 2016 Bruno S <bruno@oac.unc.edu.ar>
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
program ej9b
use precision, pr => dp

implicit none
integer                 :: i, j, k, n
real(pr), parameter     :: e = 5._pr, h = 1e-4
real(pr)                :: alpha, t0, tf, t, et, d
real(pr), dimension(:)  :: xi(4), xi_1(4), k1(4), k2(4), k3(4), k4(4), x0(4)

open(unit=10, file='ej9b_e5.dat', status='UNKNOWN', action='WRITE')
write(10, *) "#  t   q1   q2   p1   p2   e"

! x0 = (/ 2._pr, 0._pr, 0._pr, sqrt(2._pr*e - 4._pr) /)
alpha = 0.05_pr
t0 = 0._pr
tf = 15000._pr
n = int((tf-t0)/h)
x0(2) = 0  ! elijo q_2 == 0
t = 0._pr

d = 2._pr * sqrt(e)/10._pr  ! ancho barrido de espacio de fase
do k = 0, 10
    x0(3) = -sqrt(e) + real(k, pr) * d
    do j = 0, 10
        x0(1) = -sqrt(e) + real(j, pr) * d

        x0(4) = 2._pr * e - x0(3)*x0(3) - x0(1)*x0(1)
        if (x0(4) < 0) cycle
        x0(4) = sqrt(x0(4))
        xi = x0
        do i =1, n
            xi_1 = xi
            k1 = h * f(xi_1, alpha)
            k2 = h * f(xi_1 + k1*0.5_pr, alpha)
            k3 = h * f(xi_1 + k2*0.5_pr, alpha)
            k4 = h * f(xi_1 + k3, alpha)

            xi = xi_1 + (k1 + 2._pr * k2 + 2._pr * k3 + k4) / 6._pr
            t = t0 + i * h

            et = 0.5_pr*(xi(1)*xi(1)+xi(2)*xi(2)+xi(3)*xi(3)+xi(4)*xi(4)) + &
                & alpha * xi(1)*xi(1)*xi(2)*xi(2)

            if (xi(2)*xi_1(2) <= 0._pr) then
                if (xi(4) > 0._pr) then
                    write(10, *) t, xi(1), xi(2), xi(3), xi(4), et
                end if
            end if

        end do
    end do
end do
contains

function f(x, alpha)
integer(kind=8), parameter :: n = 4
real (pr) :: f(n)
real (pr), intent(in) :: x(n), alpha

    f(1) = x(3)
    f(2) = x(4)
    f(3) = -x(1) - alpha*x(1)*x(2)*x(2)
    f(4) = -x(2) - alpha*x(1)*x(1)*x(2)

end function f

end program ej9b

max_stars_repo_name: BrunoSanchez/fisicaComp

max_issues_repo_name: BrunoSanchez/fisicaComp

max_forks_repo_name: BrunoSanchez/fisicaComp


text:
c Author: Anand Patil, anand.prabhakar.patil@gmail.com

      SUBROUTINE checksymm(X,n,cs)
c Checks symmetry
cf2py intent(hide) n
cf2py intent(out) cs
      DOUBLE PRECISION X(n,n)
      INTEGER n,i,j
      LOGICAL cs

      cs = .FALSE.
      do j=1,n-1
          do i=j+1,n
              if (X(i,j).NE.X(j,i)) then
                  cs = .TRUE.
                  return
              end if
          end do
      end do

      return
      END

      SUBROUTINE chol_mvnorm(x, mu, sig, n, like, info)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(in) :: sig
cf2py double precision intent(out) :: like
cf2py integer intent(hide) :: info
cf2py threadsafe

      DOUBLE PRECISION sig(n,n), x(n), mu(n), like
      INTEGER n, info, i
      DOUBLE PRECISION twopi_N, log_detC, gd
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system
      EXTERNAL DAXPY
! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
! NB DDOT from ATLAS, compiled with gfortran 4.2 on Ubuntu Gutsy,
! was producing bad output- hence the manual dot product.

!     x <- (x-mu)
      call DAXPY(n, -1.0D0, mu, 1, x, 1)

!       mu <- x
      call DCOPY(n,x,1,mu,1)

!     x <- sig ^-1 * x
      call DPOTRS('L',n,1,sig,n,x,n,info)

      gd=0.0D0
      do i=1,n
          gd=gd+x(i)*mu(i)
      end do

!     like <- .5 dot(x,mu) (.5 (x-mu) C^{-1} (x-mu)^T)
      like = -0.5D0 * gd
!       print *, like

      twopi_N = 0.5D0 * N * dlog(2.0D0*PI)
!       print *, twopi_N

      log_detC = 0.0D0
      do i=1,n
        log_detC = log_detC + log(sig(i,i))
      enddo
!       print *, log_detC

      like = like - twopi_N - log_detC

      return
      END

      SUBROUTINE cov_mvnorm(x, mu, C, n, like, info)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(copy) :: C
cf2py double precision intent(out) :: like
cf2py integer intent(hide) :: info
cf2py threadsafe

      DOUBLE PRECISION C(n,n), x(n), mu(n), like
      INTEGER n, info
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization

!     C <- cholesky(C)
      call DPOTRF( 'L', n, C, n, info )
!       print *, C

!      If cholesky failed, puke.
       if (info .GT. 0) then
         like = -infinity
         RETURN
       endif

!     Call to chol_mvnorm
      call chol_mvnorm(x,mu,C,n,like,info)

      return
      END

      SUBROUTINE prec_mvnorm(x, mu, tau, n, like)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(copy) :: tau
cf2py double precision intent(out) :: like
cf2py threadsafe

      DOUBLE PRECISION tau(n,n), x(n), mu(n), like, gd
      INTEGER n, info, i
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)
      DOUBLE PRECISION twopi_N, log_dettau

      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DSYMV
! DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) Symmetric matrix-vector multiply
      EXTERNAL DAXPY
! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y

      twopi_N = 0.5D0 * N * dlog(2.0D0*PI)

!     x <- (x-mu)
      call DAXPY(n, -1.0D0, mu, 1, x, 1)

!       mu <- x
      call DCOPY(n,x,1,mu,1)

!     mu <- tau * x
      call DSYMV('L',n,1.0D0,tau,n,x,1,0.0D0,mu,1)

      gd = 0.0D0
      do i=1,n
          gd=gd+x(i)*mu(i)
      end do

!     like <- -.5 dot(x,mu) (.5 (x-mu) C^{-1} (x-mu)^T)
      like = -0.5D0 * gd

!      Cholesky factorize tau for the determinant.
       call DPOTRF( 'L', n, tau, n, info )

!      If cholesky failed, puke.
       if (info .GT. 0) then
         like = -infinity
         RETURN
       endif

!      Otherwise read off determinant.
       log_dettau = 0.0D0
       do i=1,n
         log_dettau = log_dettau + dlog(tau(i,i))
       enddo

       like = like - twopi_N + log_dettau

      return
      END

      SUBROUTINE blas_wishart(X,k,n,T,like)

c Wishart log-likelihood function.

cf2py double precision dimension(k,k),intent(copy) :: X,T
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,n
      DOUBLE PRECISION X(k,k),T(k,k),bx(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DSYMM
! DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) alpha*A*B + beta*C when side='l'
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization

c Check X for symmetry
      CALL checksymm(X,k,csx)
      if (csx) then
          like = -infinity
          return
      end if

c trace of T*X
!     bx <- T * X
      call DSYMM('L','L',k,k,1.0D0,T,k,x,k,0.0D0,bx,k)

c Cholesky factor T, puke if not pos def.
      call DPOTRF( 'L', k, T, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Cholesky factor X, puke if not pos def.
      call DPOTRF( 'L', k, X, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Get the trace and log-sqrt-determinants
      tbx = 0.0D0
      dx = 0.0D0
      db = 0.0D0

      do i=1,k
        tbx = tbx + bx(i,i)
        dx = dx + dlog(X(i,i))
        db = db + dlog(T(i,i))
      enddo

      if (k .GT. n) then
        like = -infinity
        RETURN
      endif

      like = 0.0D0
      like = (n - k - 1) * dx
      like = like + n * db
      like = like - 0.5D0 * tbx
      like = like - (n*k/2.0d0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)
!
      return
      END

c

      SUBROUTINE blas_wishart_cov(X,k,n,V,like)

c Wishart log-likelihood function.
c Doesn't vectorize the determinants, just the matrix multiplication.

cf2py double precision dimension(k,k),intent(copy) :: V,X
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,info, n
      DOUBLE PRECISION X(k,k),V(k,k),bx(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)
c
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system

c determinants
c Check X for symmetry
      CALL checksymm(X,k,csx)
      if (csx) then
          like = -infinity
          return
      end if

c Cholesky factorize sigma, puke if not pos def
!     V <- cholesky(V)
      call DPOTRF( 'L', k, V, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c trace of V^{-1}*X
!     bx <- X
      call DCOPY(k * k,X,1,bx,1)
!     bx <- V^{-1} * bx
      call DPOTRS('L',k,k,V,k,bx,k,info)

!     X <- cholesky(X)
      call DPOTRF( 'L', k, X, k, info )

! sqrt-log-determinant of sigma and X, and trace
      db=0.0D0
      dx=0.0D0
      tbx = 0.0D0
      do i=1,k
        db = db + dlog(V(i,i))
        dx = dx + dlog(X(i,i))
        tbx = tbx + bx(i,i)
      enddo

      if (k .GT. n) then
        like = -infinity
        RETURN
      endif

      like = (n - k - 1) * dx
      like = like - n * db
      like = like - 0.5D0*tbx
      like = like - (n*k/2.0D0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)

      return
      END

      SUBROUTINE blas_inv_wishart(X,k,n,T,like)

c Inverse Wishart log-likelihood function.

cf2py double precision dimension(k,k),intent(copy) :: X,T
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,n
      DOUBLE PRECISION X(k,k),T(k,k),TX(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DTRMM
! DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B := alpha*B*op( A )

c Check X and T for symmetry
C       CALL checksymm(X,k,csx)
C       if (csx) then
C           print *,"checksymm failed"
C           like = -infinity
C           return
C       end if

c trace of T*X^{-1}
!     TX <- T
      call DCOPY(k*k,T,1,TX,1)
!     TX <- T * X^{-1}
      call DTRMM('R','L','T','N',k,k,1.0D0,X,k,TX,k)

c Cholesky factor T, puke if not pos def.
      call DPOTRF( 'L', k, T, k, info )
      if (info .GT. 0) then
        print *,"cholesky of T failed"
        like = -infinity
        RETURN
      endif

c Cholesky factor X, puke if not pos def.
      call DPOTRF( 'L', k, X, k, info )
      if (info .GT. 0) then
        print *,"cholesky of X failed"
        like = -infinity
        RETURN
      endif

c Get the trace and log-sqrt-determinants
      tbx = 0.0D0
      dx = 0.0D0
      db = 0.0D0

      do i=1,k
        tbx = tbx + TX(i,i)
        dx = dx + dlog(X(i,i))
        db = db + dlog(T(i,i))
      enddo

      if (k .GT. n) then
        print *, "k > n"
        like = -infinity
        RETURN
      endif

      like = -1.0D0*(n + k + 1) * dx
      like = like + 0.5D0 * n * db
      like = like - 0.5D0 * tbx
      like = like - (n*k/2.0d0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)
!
      return
      END

max_stars_repo_name: CamDavidsonPilon/pymc

max_issues_repo_name: kforeman/pymc

max_forks_repo_name: kforeman/pymc


text:
C
C	$Id$
C
      DOUBLE PRECISION FUNCTION CSSIG0 (N1,N2,N,X,Y,Z,H,LIST,
     .                                LPTR,LEND,GRAD,IFLGB,
     .                                HBND,TOL,IFLGS, SIGMA,
     .                                IER)
      INTEGER N1, N2, N, LIST(*), LPTR(*), LEND(N), IFLGB,
     .        IFLGS, IER
      DOUBLE PRECISION X(N), Y(N), Z(N), H(N), GRAD(3,N),
     .                 HBND, TOL, SIGMA(*)
C
C***********************************************************
C
C                                              From SSRFPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   11/21/96
C
C   Given a triangulation of a set of nodes on the unit
C sphere, along with data values H and gradients GRAD at the
C nodes, this function determines the smallest tension fac-
C tor CSSIG0 such that the Hermite interpolatory tension
C spline H(A), defined by CSSIG0 and the endpoint values and
C directional derivatives associated with an arc N1-N2, is
C bounded (either above or below) by HBND for all A in
C (A1,A2), where (A1,A2) denotes an interval corresponding
C to the arc and A is the arc-length.
C
C On input:
C
C       N1,N2 = Nodal indexes of the endpoints of an arc for
C               which the tension factor is to be computed.
C               The indexes must be distinct and lie in the
C               range 1 to N, and if IFLGS .GE. 1, they must
C               correspond to adjacent nodes in the triangu-
C               lation.
C
C       N = Number of nodes in the triangulation.  N .GE. 3.
C
C       X,Y,Z = Arrays of length N containing coordinates of
C               the nodes.  X(I)**2 + Y(I)**2 + Z(I)**2 = 1.
C
C       H = Array of length N containing data values at the
C           nodes.  H(I) is associated with (X(I),Y(I),Z(I))
C           for I = 1 to N.
C
C       LIST,LPTR,LEND = Data structure defining the trian-
C                        gulation.  Refer to STRIPACK
C                        Subroutine CSTRMESH.
C
C       GRAD = Array dimensioned 3 by N whose columns con-
C              tain gradients at the nodes.  GRAD( ,J) must
C              be orthogonal to node J:  GRAD(1,J)*X(J) +
C              GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0.  Refer
C              to Subroutines CSGRADG, CSGRADL, and CSSMSURF.
C
C       IFLGB = Bound option indicator:
C               IFLGB = -1 if HBND is a lower bound on H.
C               IFLGB = 1 if HBND is an upper bound on H.
C
C       HBND = Bound on H.  HBND .LE. min(H1,H2) if IFLGB =
C              -1 and HBND .GE. max(H1,H2) if IFLGB = 1,
C              where H1 and H2 are the data values at the
C              endpoints of the arc N1-N2.
C
C       TOL = Tolerance whose magnitude determines how close
C             CSSIG0 is to its optimal value when nonzero
C             finite tension is necessary and sufficient to
C             satisfy the constraint.  For a lower bound,
C             CSSIG0 is chosen so that HBND .LE. HMIN .LE.
C             HBND + abs(TOL), where HMIN is the minimum
C             value of H on the arc, and for an upper bound,
C             the maximum of H satisfies HBND - abs(TOL)
C             .LE. HMAX .LE. HBND.  Thus, the constraint is
C             satisfied but possibly with more tension than
C             necessary.
C
C       IFLGS = Tension array option indicator:
C               IFLGS .LE. 0 if SIGMA is not to be used.
C               IFLGS .GE. 1 if SIGMA is to be updated by
C                            storing CSSIG0 in the appropriate
C                            locations.
C
C The above parameters are not altered by this function.
C
C       SIGMA = Dummy parameter (IFLGS .LE. 0) or array con-
C               taining tension factors associated with arcs
C               in one-to-one correspondence with LIST
C               entries (IFLGS .GE. 1).  Refer to Subroutine
C               CSGETSIG.
C
C On output:
C
C       SIGMA = Tension factor array updated with the new
C               value if and only if IFLGS .GE. 1 and IER
C               .GE. 0.
C
C       IER = Error indicator:
C             IER = 0 if no errors were encountered and the
C                     constraint can be satisfied with fin-
C                     ite tension.
C             IER = 1 if no errors were encountered but in-
C                     finite tension is required to satisfy
C                     the constraint (e.g., IFLGB = -1, HBND
C                     = H(A1), and the directional deriva-
C                     tive of H at A1 is negative).
C             IER = -1 if N1, N2, N, or IFLGB is outside its
C                      valid range.
C             IER = -2 if nodes N1 and N2 coincide or IFLGS
C                      .GE. 1 and the nodes are not adja-
C                      cent.
C             IER = -3 if HBND is outside its valid range.
C
C       CSSIG0 = Minimum tension factor defined above unless
C              IER < 0, in which case CSSIG0 = -1.  If IER
C              = 1, CSSIG0 is set to 85, resulting in an
C              approximation to the linear interpolant of
C              the endpoint values.
C
C STRIPACK module required by CSSIG0:  CSSTORE
C
C SSRFPACK modules required by CSSIG0:  CSARCLEN, CSSNHCSH
C
C Intrinsic functions called by CSSIG0:  ABS, EXP, LOG, MAX,
C                                        MIN, DBLE, SIGN,
C                                        SQRT
C
C***********************************************************
C
      DOUBLE PRECISION CSARCLEN, CSSTORE
      INTEGER LP1, LP2, LPL, LUN, NIT
      DOUBLE PRECISION A, A0, AA, AL, B, B0, BND, C, C1, C2,
     .                 COSHM, COSHMM, D, D0, D1PD2, D2,
     .                 DMAX, DSIG, E, EMS, F, F0, FMAX,
     .                 FNEG, FTOL, H1, H2, P1(3), P2(3), R,
     .                 RF, RSIG, RTOL, S, S1, S2, SBIG, SCM,
     .                 SIG, SINHM, SNEG, SSINH, SSM, STOL,
     .                 T, T0, T1, T2, TM, UN(3), UNORM
C
      DATA SBIG/85./,  LUN/-1/
      RF = DBLE(IFLGB)
      BND = HBND
C
C Print a heading.
C
      IF (LUN .GE. 0  .AND.  RF .LT. 0.) WRITE (LUN,100) N1,
     .                                   N2, BND
      IF (LUN .GE. 0  .AND.  RF .GT. 0.) WRITE (LUN,110) N1,
     .                                   N2, BND
  100 FORMAT (//1X,'CSSIG0 -- N1 =',I4,', N2 =',I4,
     .        ', LOWER BOUND = ',E15.8)
  110 FORMAT (//1X,'CSSIG0 -- N1 =',I4,', N2 =',I4,
     .        ', UPPER BOUND = ',E15.8)
C
C Test for errors and store local parameters.
C
      IER = -1
      IF (MIN(N1,N2) .LT. 1  .OR.  N1 .EQ. N2  .OR.
     .    MAX(N1,N2,3) .GT. N  .OR.  ABS(RF) .NE. 1.)
     .   GO TO 11
      IER = -2
      IF (IFLGS .GT. 0) THEN
C
C   Set LP1 and LP2 to the pointers to N2 as a neighbor of
C     N1 and N1 as a neighbor of N2, respectively.
C
        LPL = LEND(N1)
        LP1 = LPTR(LPL)
    1   IF (LIST(LP1) .EQ. N2) GO TO 2
          LP1 = LPTR(LP1)
          IF (LP1 .NE. LPL) GO TO 1
        IF (ABS(LIST(LP1)) .NE. N2) GO TO 11
C
    2   LPL = LEND(N2)
        LP2 = LPTR(LPL)
    3   IF (LIST(LP2) .EQ. N1) GO TO 4
          LP2 = LPTR(LP2)
          IF (LP2 .NE. LPL) GO TO 3
        IF (ABS(LIST(LP2)) .NE. N1) GO TO 11
      ENDIF
C
C Store nodal coordinates P1 and P2, compute arc-length AL
C   and unit normal UN = (P1 X P2)/UNORM, and test for
C   coincident nodes.
C
    4 P1(1) = X(N1)
      P1(2) = Y(N1)
      P1(3) = Z(N1)
      P2(1) = X(N2)
      P2(2) = Y(N2)
      P2(3) = Z(N2)
      AL = CSARCLEN(P1,P2)
      UN(1) = P1(2)*P2(3) - P1(3)*P2(2)
      UN(2) = P1(3)*P2(1) - P1(1)*P2(3)
      UN(3) = P1(1)*P2(2) - P1(2)*P2(1)
      UNORM = SQRT(UN(1)*UN(1) + UN(2)*UN(2) + UN(3)*UN(3))
      IF (UNORM .EQ. 0.  .OR.  AL .EQ. 0.) GO TO 11
C
C Store endpoint data values and test for valid constraint.
C
      H1 = H(N1)
      H2 = H(N2)
      IER = -3
      IF ((RF .LT. 0.  .AND.  MIN(H1,H2) .LT. BND)  .OR.
     .    (RF .GT. 0.  .AND.  BND .LT. MAX(H1,H2)))
     .   GO TO 11
C
C Compute scaled directional derivatives S1,S2 at the end-
C   points (for the direction N1->N2) and test for infinite
C   tension required.
C
      S1 = AL*(GRAD(1,N1)*P2(1) + GRAD(2,N1)*P2(2) +
     .         GRAD(3,N1)*P2(3))/UNORM
      S2 = -AL*(GRAD(1,N2)*P1(1) + GRAD(2,N2)*P1(2) +
     .          GRAD(3,N2)*P1(3))/UNORM
      IER = 1
      SIG = SBIG
      IF ((H1 .EQ. BND  .AND.  RF*S1 .GT. 0.)  .OR.
     .    (H2 .EQ. BND  .AND.  RF*S2 .LT. 0.)) GO TO 10
C
C Test for SIG = 0 sufficient.
C
      IER = 0
      SIG = 0.
      IF (RF*S1 .LE. 0.  .AND.  RF*S2 .GE. 0.) GO TO 10
C
C   Compute first difference S and coefficients A0 and B0
C     of the Hermite cubic interpolant H0(A) = H2 - (S2*R +
C     B0*R**2 + (A0/3)*R**3), where R(A) = (A2-A)/AL.
C
      S = H2 - H1
      T0 = 3.*S - S1 - S2
      A0 = 3.*(S-T0)
      B0 = T0 - S2
      D0 = T0*T0 - S1*S2
C
C   H0 has local extrema in (A1,A2) iff S1*S2 < 0 or
C     (T0*(S1+S2) < 0 and D0 .GE. 0).
C
      IF (S1*S2 .GE. 0.  .AND.  (T0*(S1+S2) .GE. 0.  .OR.
     .    D0 .LT. 0.)) GO TO 10
      IF (A0 .EQ. 0.) THEN
C
C   H0 is quadratic and has an extremum at R = -S2/(2*B0).
C     H0(R) = H2 + S2**2/(4*B0).  Note that A0 = 0 implies
C     2*B0 = S1-S2, and S1*S2 < 0 implies B0 .NE. 0.
C     Also, the extremum is a min iff HBND is a lower bound.
C
        F0 = (BND - H2 - S2*S2/(4.*B0))*RF
      ELSE
C
C   A0 .NE. 0 and H0 has extrema at R = (-B0 +/- SQRT(D0))/
C     A0 = S2/(-B0 -/+ SQRT(D0)), where the negative root
C     corresponds to a min.  The expression for R is chosen
C     to avoid cancellation error.  H0(R) = H2 + (S2*B0 +
C     2*D0*R)/(3*A0).
C
        T = -B0 - SIGN(SQRT(D0),B0)
        R = T/A0
        IF (RF*B0 .GT. 0.) R = S2/T
        F0 = (BND - H2 - (S2*B0+2.*D0*R)/(3.*A0))*RF
      ENDIF
C
C   F0 .GE. 0 iff SIG = 0 is sufficient to satisfy the
C     constraint.
C
      IF (F0 .GE. 0.) GO TO 10
C
C Find a zero of F(SIG) = (BND-H(R))*RF where the derivative
C   of H, HP, vanishes at R.  F is a nondecreasing function,
C   F(0) < 0, and F = FMAX for SIG sufficiently large.
C
C Initialize parameters for the secant method.  The method
C   uses three points:  (SG0,F0), (SIG,F), and (SNEG,FNEG),
C   where SG0 and SNEG are defined implicitly by DSIG = SIG
C   - SG0 and DMAX = SIG - SNEG.  SG0 is initially zero and
C   SNEG is initialized to a sufficiently large value that
C   FNEG > 0.  This value is used only if the initial value
C   of F is negative.
C
      FMAX = MAX(1.D-3,MIN(ABS(H1-BND),ABS(H2-BND)))
      T = MAX(ABS(H1-BND),ABS(H2-BND))
      SIG = MAX(ABS(S1),ABS(S2))/T
      DMAX = SIG*(1.-T/FMAX)
      SNEG = SIG - DMAX
      IF (LUN .GE. 0) WRITE (LUN,120) SIG, SNEG, F0, FMAX
  120 FORMAT (1X,8X,'SIG = ',E15.8,', SNEG = ',E15.8/
     .        1X,9X,'F0 = ',E15.8,', FMAX = ',E15.8/)
      DSIG = SIG
      FNEG = FMAX
      D2 = S2 - S
      D1PD2 = S2 - S1
      NIT = 0
C
C Compute an absolute tolerance FTOL = abs(TOL) and a
C   relative tolerance RTOL = 100*Macheps.
C
      FTOL = ABS(TOL)
      RTOL = 1.
    5 RTOL = RTOL/2.
        IF (CSSTORE(RTOL+1.) .GT. 1.) GO TO 5
      RTOL = RTOL*200.
C
C Top of loop:  compute F.
C
    6 EMS = EXP(-SIG)
      IF (SIG .LE. .5) THEN
C
C   Use approximations designed to avoid cancellation error
C     (associated with small SIG) in the modified hyperbolic
C     functions.
C
        CALL CSSNHCSH (SIG, SINHM,COSHM,COSHMM)
        C1 = SIG*COSHM*D2 - SINHM*D1PD2
        C2 = SIG*(SINHM+SIG)*D2 - COSHM*D1PD2
        A = C2 - C1
        AA = A/EMS
        E = SIG*SINHM - COSHMM - COSHMM
      ELSE
C
C   Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid
C     overflow.
C
        TM = 1. - EMS
        SSINH = TM*(1.+EMS)
        SSM = SSINH - 2.*SIG*EMS
        SCM = TM*TM
        C1 = SIG*SCM*D2 - SSM*D1PD2
        C2 = SIG*SSINH*D2 - SCM*D1PD2
        AA = 2.*(SIG*TM*D2 + (TM-SIG)*D1PD2)
        A = EMS*AA
        E = SIG*SSINH - SCM - SCM
      ENDIF
C
C   HP(R) = (S2 - (C1*sinh(SIG*R) - C2*coshm(SIG*R))/E)/DT
C     = 0 for ESR = (-B +/- sqrt(D))/A = C/(-B -/+ sqrt(D))
C     where ESR = exp(SIG*R), A = C2-C1, D = B**2 - A*C, and
C     B and C are defined below.
C
      B = E*S2 - C2
      C = C2 + C1
      D = B*B - A*C
      F = 0.
      IF (AA*C .EQ. 0.  .AND.  B .EQ. 0.) GO TO 7
      F = FMAX
      IF (D .LT. 0.) GO TO 7
      T1 = SQRT(D)
      T = -B - SIGN(T1,B)
      RSIG = 0.
      IF (RF*B .LT. 0.  .AND.  AA .NE. 0.) THEN
        IF (T/AA .GT. 0.) RSIG = SIG + LOG(T/AA)
      ENDIF
      IF ((RF*B .GT. 0.  .OR.  AA .EQ. 0.)  .AND.
     .    C/T .GT. 0.) RSIG = LOG(C/T)
      IF ((RSIG .LE. 0.  .OR.  RSIG .GE. SIG)  .AND.
     .    B .NE. 0.) GO TO 7
C
C   H(R) = H2 - (B*SIG*R + C1 + RF*sqrt(D))/(SIG*E).
C
      F = (BND - H2 + (B*RSIG+C1+RF*T1)/(SIG*E))*RF
C
C   Update the number of iterations NIT.
C
    7 NIT = NIT + 1
      IF (LUN .GE. 0) WRITE (LUN,130) NIT, SIG, F
  130 FORMAT (1X,3X,I2,' -- SIG = ',E15.8,', F = ',
     .        E15.8)
      IF (F0*F .LT. 0.) THEN
C
C   F0*F < 0.  Update (SNEG,FNEG) to (SG0,F0) so that F and
C     FNEG always have opposite signs.  If SIG is closer to
C     SNEG than SG0, then swap (SNEG,FNEG) with (SG0,F0).
C
        T1 = DMAX
        T2 = FNEG
        DMAX = DSIG
        FNEG = F0
        IF (ABS(DSIG) .GT. ABS(T1)) THEN
C
          DSIG = T1
          F0 = T2
        ENDIF
      ENDIF
C
C   Test for convergence.
C
      STOL = RTOL*SIG
      IF (ABS(DMAX) .LE. STOL  .OR.  (F .GE. 0.  .AND.
     .    F .LE. FTOL)  .OR.  ABS(F) .LE. RTOL) GO TO 10
C
C   Test for F0 = F = FMAX or F < 0 on the first iteration.
C
      IF (F0 .NE. F  .AND.  (NIT .GT. 1  .OR.  F .GT. 0.))
     .   GO TO 9
C
C   F*F0 > 0 and either the new estimate would be outside
C     of the bracketing interval of length abs(DMAX) or
C     F < 0 on the first iteration.  Reset (SG0,F0) to
C     (SNEG,FNEG).
C
    8 DSIG = DMAX
      F0 = FNEG
C
C   Compute the change in SIG by linear interpolation
C     between (SG0,F0) and (SIG,F).
C
    9 DSIG = -F*DSIG/(F-F0)
      IF (LUN .GE. 0) WRITE (LUN,140) DSIG
  140 FORMAT (1X,8X,'DSIG = ',E15.8)
      IF ( ABS(DSIG) .GT. ABS(DMAX)  .OR.
     .     DSIG*DMAX .GT. 0. ) GO TO 8
C
C   Restrict the step-size such that abs(DSIG) .GE. STOL/2.
C     Note that DSIG and DMAX have opposite signs.
C
      IF (ABS(DSIG) .LT. STOL/2.) DSIG = -SIGN(STOL/2.,DMAX)
C
C   Bottom of loop:  Update SIG, DMAX, and F0.
C
      SIG = SIG + DSIG
      DMAX = DMAX + DSIG
      F0 = F
      GO TO 6
C
C No errors encountered.
C
   10 CSSIG0 = SIG
      IF (IFLGS .LE. 0) RETURN
      SIGMA(LP1) = SIG
      SIGMA(LP2) = SIG
      RETURN
C
C Error termination.
C
   11 CSSIG0 = -1.
      RETURN
      END

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text: Hey, where can I find some of these guys?  I met you at the colleges, and would love to get in touch w your dads friend, even if I have to wait a few months...   ErikMadden


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE INBRST (MS, MR, N5, N6, N21, N23, JJ, IMTRL, JC, IIN,
     &   IFOUND, IBARST, JMAT, JCENT, NLPB, JFLINE, JLLIST, LINKB,
     &   LINKM, NHOLDM, IHOLDM, NHOLDB, IHOLDB, MERGE, NOROOM)
C***********************************************************************

C  SUBROUTINE INBRST = INPUTS A BAR SET INTO THE DATABASE

C***********************************************************************

      DIMENSION IBARST (MS), JMAT (MS), JCENT (MS), NLPB (MS)
      DIMENSION JFLINE (MS)
      DIMENSION JLLIST (3 * MS), LINKB (2, MS), LINKM (2, MS + MR)
      DIMENSION IHOLDM (2,  (MS + MR)), IHOLDB (2, MS)
      DIMENSION IIN (IFOUND)

      LOGICAL MERGE, NOROOM, ADDLNK

      IZ = 0
      NOROOM = .TRUE.
      N22 = 0

C  ZERO OUT THE LINK ARRAY IF NEEDED

      IF (JJ .GT. N21) THEN
         N21 = JJ

C  FIND THE CORRECT BAR SET NUMBER IF MERGING

C  SET UP POINTERS FOR MERGING DATA

      ELSEIF (MERGE) THEN
         JHOLD = JJ
         CALL LTSORT (MS, LINKB, JJ, IPNTR, ADDLNK)
         IF (IPNTR .GT. 0) THEN
            IF (JHOLD .GT. NHOLDB)NHOLDB = JHOLD
            CALL LTSORT (MS, IHOLDB, JHOLD, IPNTR, ADDLNK)
            IF (IPNTR .GT. 0) THEN
               JJ = IPNTR
            ELSE
               JJ = N22 + 1
               N22 = JJ
               WRITE ( * , 10000)JHOLD, JJ
               ADDLNK = .TRUE.
               CALL LTSORT (MS, IHOLDB, JHOLD, JJ, ADDLNK)
            ENDIF
         ENDIF
      ENDIF

C  INPUT THE BAR SET DATA INTO THE DATABASE

      N5 = N5 + 1
      J = N5
      IF (J .GT. MS)RETURN
      ADDLNK = .TRUE.
      CALL LTSORT (MS, LINKB, JJ, J, ADDLNK)
      IBARST (J) = JJ
      JCENT (J) = JC
      JFLINE (J) = N6 + 1
      DO 100 I = 1, IFOUND
         JJ = IIN (I)
         IF (JJ .EQ. 0)GOTO 110
         N6 = N6 + 1
         IF (N6 .GT. MS * 3)RETURN
         JLLIST (N6) = JJ
  100 CONTINUE
  110 CONTINUE
      NLPB (J) = N6 - JFLINE (J) + 1
      IF (NLPB (J) .LT. 1) THEN
         WRITE ( * , 10010)J
         CALL LTSORT (MS, LINKB, JJ, IZ, ADDLNK)
      ENDIF
      ADDLNK = .FALSE.

C  LINK UP THE MATERIAL

C  ZERO THE LINK ARRAY IF NEEDED

      IF (IMTRL .GT. N23) THEN
         N23 = IMTRL

C  SET UP POINTERS FOR MERGING DATA

      ELSEIF (MERGE) THEN
         JHOLD = IMTRL
         CALL LTSORT (MS + MR, LINKM, IMTRL, IPNTR, ADDLNK)
         IF (IPNTR .NE. 0) THEN
            IF (JHOLD .GT. NHOLDM)NHOLDM = JHOLD
            ADDLNK = .FALSE.
            CALL LTSORT ( (MS + MR), IHOLDM, JHOLD, IPNTR, ADDLNK)
            IF (IPNTR .GT. 0) THEN
               IMTRL = IPNTR
            ELSE
               IMTRL = N23 + 1
               N23 = IMTRL
               WRITE ( * , 10010)JHOLD, IMTRL
               ADDLNK = .TRUE.
               CALL LTSORT ( (MS + MR), IHOLDM, JHOLD, IMTRL, ADDLNK)
            ENDIF
         ENDIF
      ENDIF

C  ADD THE MATERIAL INTO THE DATABASE

      NOROOM = .FALSE.
      ADDLNK = .FALSE.
      CALL LTSORT (MS + MR, LINKM, IMTRL, IPNTR, ADDLNK)
      ADDLNK = .TRUE.
      IF (IPNTR .GT. 0) THEN
         CALL MESAGE (' ')
         WRITE ( * , 10020)IMTRL, IBARST (J)
         CALL LTSORT (MS, LINKB, IBARST (J), IZ, ADDLNK)
         RETURN
      ELSEIF (IPNTR .EQ. 0) THEN
         IMINUS =  - 1
         CALL LTSORT (MS + MR, LINKM, IMTRL, IMINUS, ADDLNK)
      ENDIF
      JMAT (J) = IMTRL

      RETURN

10000 FORMAT ('   OLD BAR SET NO:', I5, '   TO NEW BAR SET NO:', I5)
10010 FORMAT (' BAR SET:', I5, ' HAS LESS THAN ONE LINE',  / ,
     &   ' THIS BAR SET WILL NOT BE INPUT INTO DATABASE')
10020 FORMAT (' MATERIAL:', I5, ' FOR BAR SET:', I5,
     &   ' HAS BEEN DESIGNATED',
     &   / , ' AS A REGION  (4 NODE ELEMENT) MATERIAL.',  / ,
     &   ' ELEMENTS WITH 2 AND 4 NODES CANNOT SHARE MATERIAL ID''S',
     &   / , ' THIS BAR SET WILL NOT BE INPUT INTO DATABASE')
      END

max_stars_repo_name: utkarshayachit/seacas

max_issues_repo_name: utkarshayachit/seacas

max_forks_repo_name: utkarshayachit/seacas