program main
       
       integer nd, nx
       parameter (nd=8, nx=32)

       double precision margin
       parameter (margin = 1.)

       double precision begintime, endtime
       parameter (begintime = 0., endtime = 1.)

       integer dumps
       parameter (dumps = 10)

       integer    ll, mm
       parameter  (ll=2, mm=2)

       double precision time
       double precision dd, dx, du
       integer    i, it, nt, scoop

       double precision x(0:nx)
       double precision q(-(nd+1):(nd+1))
       double precision p(-(nd+1):(nd+1))
       double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
       double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
       integer     mask(-(nd+1):(nd+1),-(nd+1):(nd+1))

       integer aa, bb, cc, dd1

       INTERFACE
         subroutine interpolate(gn, q, p, nx, nd, dd, i)
       
       double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
       double precision q(-(nd+1):(nd+1))
       double precision p(-(nd+1):(nd+1))
       integer nx, nd, i
       double precision dd

       end subroutine interpolate
       
       subroutine gmiddle(go,gn,x,q,p,nd,nx,dd,dx,du,i)
       double precision x(0:nx)
       double precision q(-(nd+1):(nd+1))
       double precision p(-(nd+1):(nd+1))
       double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
       double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
       integer          nd, nx
       double precision dd, dx, du
       integer          i
       end subroutine gmiddle    
       
       END INTERFACE

       open(24,file='pfvals',status='new')
       call setup(mask, x, q, p, nd, nx , dd, dx, du, ll,
     & mm, margin, time, begintime, endtime, nt
     & , dumps, scoop)

       call initial(gn, x, q, p, nd, nx, time, ll, mm)

       do i=0,nx
       call interpolate(gn, q, p, nx, nd, dd, i)
       enddo

       go=gn

       do it=1,nt
       time=time+du
       i=1
       call gstart(go,gn,x,q,p,nd,nx,dd,dx,du)

      
       call interpolate(gn,q,p,nx,nd,dd,i)
       do i=2, nx-1
       call gmiddle(go,gn,x,q,p,nd,nx,dd,dx,du,i)
       call interpolate(gn,q,p,nx,nd,dd,i)
       enddo
       i=nx
       call gend(go,gn,x,q,p,nd,nx,dd,dx,du)
       call interpolate(gn,q,p,nx,nd,dd,i)
       call norm(gn,x,q,p,nd,nx,dd,dx,time,dumps,ll,mm) 
 
c       open(22,file='gnvals2',status='new')
c       do aa=-(nd+1), (nd+1)
c       do bb=-(nd+1),(nd+1)
c       do cc=0,nx
c       do dd1=0,1
     
c           write(22,*) aa,bb,cc,dd1,gn(aa,bb,cc,dd1)
c       enddo
c       enddo
c       enddo
c       enddo

c       call flush(22)
        go=gn
        enddo
       call flush(24)


       contains


   
       

      subroutine setup(mask, x, q, p, nd, nx, dd, dx, du, ll, mm, 
     &                 margin, time, begintime, endtime, nt, 
     &                 dumps, scoop)

      integer mask(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer nd, nx
      double precision dd, dx, du
      integer ll, mm
      double precision margin, time, begintime, endtime
      integer nt, dumps, scoop

      integer i

C     Fill-in grid coordinate functions
C     ND is the real grid dimension, the grid is enlarged by one,
C     so that it holds the interpolated values as well.

      dd = 1. / dble(nd)
      dx = 1. / dble(nx)


      do i = 0, nx 
         x(i)  = i * dx

      end do 

      do i = -(nd+1), (nd+1)
         q(i)  = i * dd
         p(i)  = i * dd
      end do

C     Figure out the number of iterations (nt)
C     Calculate the timestep (du), and the sampling density (scoop).

      time = begintime
      du = margin * dx * dd ** 2

      nt = int((endtime - du) / du) + 1
      if (nt .lt. dumps) dumps = nt
      if (dumps .gt. 0) then
        scoop = int(nt / dumps)
      else
        scoop = 1
      end if

      if (dumps .gt. 0) then
        open (11,file='dat/run.log' ,status='new')
        write (11,*) 'Run Parameters'
        write (11,*) '-------------------------------------'
        write (11,*) 'Grid size            : ', nx, ' x ', nd, ' x ', nd
        write (11,*) 'Initial time         : ', time
        write (11,*) 'Final time           : ', endtime
        write (11,*) 'Number of iterations : ', nt
        write (11,*) 'Number of data dumps : ', dumps
        write (11,*) ' '
        call flush(11) 
      end if

C     Set mask array for interpolation between stereographic patches

c      call setmask(mask, q, p, nd)

      return
      end subroutine setup

        subroutine initial(gn, x, q, p, nd, nx, time, ll, mm)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      double precision time
      integer          nd, nx, ll, mm

      integer          i, k, l
      double precision u

C     Fill in initial data for all x's, both patches, all angles

      u = time

      do i = 0, nx 
        do l = -nd, nd
          do k = -nd, nd
            gn(k,l,i,0) = fg(ll,mm,u,q(k),p(l),x(i),0)
            gn(k,l,i,1) = fg(ll,mm,u,q(k),p(l),x(i),1)
         end do
       end do
      end do

      return
      end subroutine initial

      double precision function fg(ll, mm, u, q, p, x, D)
      double precision u, q, p, x, y
      double precision lambda
c      double precision cos1, ylm, cos2, phi
      double precision phi
      integer ll,mm, D, l, m

C     Choose overall factor

      lambda = 1.
      l = ll
      m = mm

C     Calculate the cosine for this point and this patch

      y = cos1(q,p,D)
      phi = acos(cos2(q,p,D))
      if(abs(cos2(q,p,D)).gt.1) then
      print*, q, p, cos2(q,p,D)
      end if

      fg = lambda * 2**(l+1) *  x ** (l+1) * ylm(l, m, y, phi)
     &   / ((u + 1.) ** (l + 1) * (u * (1 - x) + x + 1.) ** (l+1))

      return 
      end function fg

      double precision function legendre(l,y)

C     Legendre Polynomials

      integer l, lc
      double precision y
      double precision pnext, tmp0, tmp1

C     l = 0 and 1 polynomials 

      tmp0 = 1.
      tmp1 = y

      if(l.eq.0) then
         legendre = tmp0
      else if (l.eq.1) then 
         legendre = tmp1
      else 

C     Higher order recursive formula

         do lc = 2, l
            pnext = (y * (2*lc-1) * tmp1 - (lc-1) * tmp0) / lc
            tmp0 = tmp1
            tmp1 = pnext
         end do

         legendre = pnext

      endif
      return
      end function legendre

      double precision function ylm(l, m, x, phi)

C     Spherical Harmonics
C     Re(Ylm) (The real part)

      integer l, m
      double precision x, phi
      double precision pi, norm
c      double precision plgndr, factrl

      pi = 3.14159265358979323846

      norm = sqrt((2. * l + 1.) / (4. * pi) 
     &       * factrl(l - m) / factrl(l + m))

      ylm = norm * plgndr(l, m, x) * cos(m * phi)

      return
      end function ylm

      double precision function plgndr(l,m,x)
      integer l,m
      double precision x
      integer i,ll
      double precision fact, pll, pmm, pmmp1, somx2

      if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.) then
      stop 'bad arguments to plgndr'
      end if

      pmm=1.
      if(m.gt.0) then
      somx2= sqrt((1. -x)*(1. +x))
      fact=1.
      do i=1,m
          pmm=-pmm*fact*somx2
          fact=fact+2
      enddo
      endif
      
      if(l.eq.m) then
          plgndr=pmm
      else
        pmmp1=x*(2*m+1)*pmm
      if(l.eq.m+1) then
        plgndr=pmmp1
      else
        do ll=m+2,l
           pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
           pmm=pmmp1
           pmmp1=pll
        enddo
       plgndr=pll
       endif
       endif
       return 
       end  function plgndr

      double precision function factrl(n)

C     Factorial

      integer n
c      double precision gammln
      integer j,ntop
      double precision a(33)
      save ntop,a
      data ntop,a(1)/0,1./

      if (n.lt.0) then
          stop 'negative factorial'
      else if (n.le.ntop) then
          factrl=a(n+1)
      else if (n.le.32) then
          do j=ntop+1,n
            a(j+1)=j*a(j)
          end do
          ntop=n
          factrl=a(n+1)
      else
      factrl = dexp(gammln(dble(n)+1.0D0))
      endif
      return
      end function factrl

      double precision function gammln(xx)
      double precision xx
      double precision half,one,fpf
      parameter (half=0.5d0, one=1.0d0, fpf=5.5d0)
      integer j
      double precision ser,stp,tmp,x,cof(6)
      save cof,stp
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     & -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/

      x = xx - one
      tmp = x + fpf
      tmp = (x+half) * log(tmp) - tmp
      ser = one

      do j=1,6
         x=x+one
         ser=ser+cof(j)/x
      end do
      gammln = tmp+log(stp*ser)

      return
      end function gammln

      double precision function p1(u,q,p,x,D)
      double precision u, q, p, x
      double precision lambda
c      double precision cos2, sin1
      integer D
      
      lambda = 1.

      p1 = lambda * 4 * sin1(q,p,D) * cos2(q,p,D) * x**2 /
     &    ( (-x - u + u * x - 1)**2 * (u + 1)**2)


      return 
      end function p1

      double precision function sin1(q,p,D)
      double precision q, p  
      integer D
      sin1 = 2. * sqrt(q*q + p*p) / (1. + q*q + p*p) 
      if (sin1.gt.1) print*, 'Error in sin1'
      end  function sin1

      double precision function sin2(q,p,D)
      double precision q, p  
      integer D
      if(q.eq.0.and.p.eq.0) then
        sin2 = 0.
      else
          if(D.eq.0) then 
               sin2 =   p / sqrt(q*q + p*p)
          else if(D.eq.1) then
               sin2 = - p / sqrt(q*q + p*p)
          end if 
      end if
      if (sin2.gt.1) print*, 'Error in sin2'
      end function sin2

      double precision function cos1(q,p,D)
      double precision q, p  
      integer D
      if(D.eq.0) then 
         cos1 = (q*q + p*p -1.) / (1. + q*q + p*p)
      else if(D.eq.1) then
         cos1 = (1. - q*q - p*p) / (1. + q*q + p*p) 
      end if
      if (cos1.gt.1) print*, 'Error in cos1'
      end function cos1

      double precision function cos2(q,p,D)
      double precision q, p  
      integer D
      if(q.eq.0.and.p.eq.0) then
        cos2 = 0.
      else if(p.eq.0.and.q.gt.0) then
        cos2 = 1.
      else if(p.eq.0.and.q.lt.0) then
        cos2 = - 1.
      else 
        cos2 = q / sqrt(q*q + p*p)
      end if
      if (cos2.gt.1) cos2 = 1.
      end function cos2
      subroutine copy(go, gn, nd, nx)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      integer          nd, nx

      integer          i, k, l

C     Copy the "new" array into the "old" one. We could use a
C     "canned" routine from a library to do this chore.

      do i = 0, nx
        do l = -(nd+1), (nd+1)
          do k = -(nd+1), (nd+1)
             go(k,l,i,0) = gn(k,l,i,0)
             go(k,l,i,1) = gn(k,l,i,1)
          end do
        end do
      end do

      return
      end subroutine copy

      subroutine gstart(go, gn, x, q, p, nd, nx, dd, dx, du)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, du

      integer          D, i, k, l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs
      double precision rp, rq, rr, rs, rc
      double precision gp, gq, gr, gs
      double precision l2, l2im1, l2i, factor, area
      double precision gmean

C     Calculate the characteristic deviation 

      i = 1
      dold = 0.5  * du * (1. - x(i-1)) ** 2 
      dnew = 0.25 * du * (1. - x(i))   ** 2 
      d1 = dold / dx 
      d2 = dnew / dx 

      xp = x(i-1) 
      xq = x(i)   - dnew
      xr = x(i-1) + dold
      xs = x(i)   + dnew

      rp = xp / (1. - xp)
      rq = xq / (1. - xq)
      rr = xr / (1. - xr)
      rs = xs / (1. - xs)

      rc = 0.5 * (rp + rs)

C     area = \int_\Sigma du dr

      area = 0.5 * du * (rq - rp + rs - rr)
 
C     Calculate gn for both patches independently

      do D = 0, 1
        gmean = 0.
        do l = -nd, nd
          do k = -nd, nd

            factor = - 0.25 * (1 + q(k) ** 2 + p(l) ** 2) ** 2
            l2im1 = 0.
            l2i = factor * (1. / dd ** 2) * 
     &            (go(k+1,l,i,D)
     &           + go(k-1,l,i,D) + go(k,l+1,i,D)
     &           + go(k,l-1,i,D) - 4. * go(k,l,i,D)) 
            l2 = 0.5 * (l2im1 + l2i)

C           Interpolate the g = r * f values

            gp = 0.

            gr = (1. - d1) * go(k,l,i-1,D) + d1 * go(k,l,i,D)

            gs = (1. - d2) * go(k,l,i,D)   + d2 * go(k,l,i+1,D)

            gq =  gp + gs - gr - 0.5 * area * l2 / rc ** 2

C           Extrapolate to obtain the new gamma grid value

            gn(k,l,i,D) = (gq - gn(k,l,i-1,D) * d2) / (1. - d2)
            gmean = gmean + gn(k,l,i,D)
          end do
        end do

        gmean = gmean / (2.* nd + 1.) ** 2

C       Average out the angular dependence

      forall (l=-nd:nd,k=-nd:nd) gn(k,l,i,D) = gmean

      end do

      return
      end subroutine gstart

      subroutine gend(go, gn, x, q, p, nd, nx, dd, dx, du)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, du

      integer          D, i, k, l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs
      double precision rp,  rr
      double precision gp, gq, gr, gs
      double precision l2, l2im1, l2i, factor, integral

C     Calculate the characteristic deviation 

      i = nx
      dold = 0.25 * du * (1. - x(i-1)) ** 2
      dnew = 0.
      d1 = dold / dx 
      d2 = dnew / dx 

      xp = x(i-1) - dold
      xq = 1.
      xr = x(i-1) + dold
      xs = 1.

      rp = xp / (1. - xp)
      rr = xr / (1. - xr)

C     integral = \int_\Sigma du dr 1/r^2

      integral = 2. * log(rr / rp)
 
C     Calculate gn for both patches independently

      do D = 0, 1
        do l = -nd, nd
          do k = -nd, nd

C           Calculate the Laplacian
       
            factor = - 0.25 * (1 + q(k) ** 2 + p(l) ** 2) ** 2

            l2im1 = factor * (1. / dd ** 2)
     &            * (gn(k+1,l,i-1,D) + gn(k-1,l,i-1,D) 
     &             + gn(k,l+1,i-1,D) + gn(k,l-1,i-1,D) 
     &             - 4. * gn(k,l,i-1,D)) 

            l2i = factor * (1. / dd ** 2)  
     &          * (go(k+1,l,i,D) + go(k-1,l,i,D) 
     &           + go(k,l+1,i,D) + go(k,l-1,i,D)
     &           - 4. * go(k,l,i,D)) 

            l2 = 0.5 * (l2im1 + l2i)

C           Interpolate the g = r * f values

            gp = (1. - d1) * gn(k,l,i-1,D) + d1 * gn(k,l,i-2,D)

            gr = (1. - d1) * go(k,l,i-1,D) + d1 * go(k,l,i,D)

            gs = go(k,l,i,D)

C           Calculate the gp box value

            gq = gp + gs - gr - 0.5 * integral * l2

C           Extrapolate to obtain the new gamma grid value

            gn(k,l,i,D) = (gq - gn(k,l,i-1,D) * d2) / (1. - d2)

          end do
        end do
      end do

      return
      end subroutine gend
      subroutine norm(gn, x, q, p, nd, nx, dd, dx, time, dumps, ll, mm)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, time
      integer          dumps, ll, mm

      integer          i, k, l, D
      double precision ex, nu, diff, l1, l2, linf
      integer          mk, ml, mi

C     Calculate the global l1, l2 and linfinity norms.

      l1 = 0.
      l2 = 0.
      linf = 0.

      mk = 0
      ml = 0
      mi = 0

      D = 0
      do i = 0, nx
        do l = -nd, nd
          do k = -nd, nd
           ex = fg(ll,mm,time,q(k),p(l),x(i),D)
           nu = gn(k,l,i,D) 
           diff = abs(ex - nu)
           l1 = l1 + diff
           l2 = l2 + diff * diff
           if (diff.ge.linf) then
               mk = k
               mi = i
               ml = l
               linf = diff
           end if
         end do 
        end do 
      end do

      l2 = sqrt(l2) / sqrt(dble(nx * nd * nd))
      l1 = l1 / (dble(nx * nd * nd))

      if (dumps .gt. 0) then
        open(14, file = 'dat/norml1.dat', status = 'unknown')
        open(15, file = 'dat/norml2.dat', status = 'unknown')
        open(16, file = 'dat/normli.dat', status = 'unknown')
        write (14,10) time, l1
        write (15,10) time, l2
        write (16,20) time, linf, mk, ml, mi
        call flush (14)
        call flush (15)
        call flush (16)
      end if

   10 format(2E23.12E3)
   20 format(2E23.12E3,3I5)

      return
      end subroutine norm  
      
      end program main