      subroutine xskdline (ccheck,ear,ne,param,ifl,photar,photer)

      implicit none

      integer ni,ngp,nc,ng,nr,np,nb
      parameter (ni=2,ngp=5,nc=2,ng=100,nr=100,np=5,nb=100)

      integer ne, ifl
      real ear(0:ne), param(9), photar(ne), photer(ne)
      logical ccheck

c This model calculates the line profile from a Keplerian Accretion disk around
c both Schwarzschild & Kerr (arbitrary spin) black holes.
c
c It requires the use of an extensive set of transfer functions, which must be
c downloaded from the XSPEC website. These transfer functions should be placed
c in $LHEASOFT/../spectral/xspec/manager/XSKDLINE/
c
c It is not called directly from XSPEC, rather, it is called by the files
c xskdlN (additive model), xskdcN (convolution model, via xskdconv). Here,
c N is either 0 or 1 and denotes different modes of operation for the model.
c For further information how each of the modes works, consult xskdlN.f,
c xskdcN.f, xskdconv.f as appropriate.
c
c Finally, the user should note that chatter>=12 will case the model to dump
c the current line parameters to the screen

      double precision rval
      double precision r1,r2
      double precision a,as,obs
      double precision risco,rms
      double precision ec,emin,emax
      double precision rmin,rmax,index
      integer rtype,mtype,ftype,fts,nes
      character*80 buff
      logical qfirst

      double precision de(0:nb-1)
      double precision e(0:nb-1),n(nb)

      logical nmask(0:nb-1)
      logical check1,check2,check4

      integer loc,i

      save nmask,n,e,de,r1,r2
      save ec,a,as,obs,rmin,rmax,index
      save rtype,mtype,ftype,fts,nes
      save qfirst

      external rms
      
      DATA qfirst/.TRUE./



      ftype = int(param(9))

      If (qfirst) then

         nes = ne

         fts = ftype

         as = dble(param(2))

         If (ftype.EQ.0) then

            emin = dble(ear(0))
            emax = dble(ear(ne-1))
            ec = dble(param(1))
            a = dble(param(2))
            rmin = dble(param(3))
            rmax = dble(param(4))
            obs = dble(param(5))
            index = dble(param(6))
            rtype = int(param(7))
            mtype = int(param(8))

         Else

            emin = dble(ear(0))
            emax = dble(ear(ne-1))
            ec = dble(param(1))
            a = dble(param(2))
            rmin = dble(param(3))
            rmax = dble(param(4))
            obs = dble(param(5))
            index = dble(param(6))
            rtype = int(param(7))
            mtype = int(param(8))

         End If

         write(*,*)'This is the XSKDLine model by Kris Beckwith'


      End If


      check1=.FALSE.

      If (qfirst) check1=.TRUE.

      If (ftype.NE.fts) check1=.TRUE.

      If (ftype.EQ.0) then

         If (rmin.NE.dble(param(3))) check1=.TRUE.
         If (rmax.NE.dble(param(4))) check1=.TRUE.
         If (obs.NE.dble(param(5))) check1=.TRUE.
         If (index.NE.dble(param(6))) check1=.TRUE.
         If (rtype.NE.int(param(7))) check1=.TRUE.
         If (mtype.NE.int(param(8))) check1=.TRUE.

      Else if ((ftype.EQ.1).OR.(ftype.EQ.2)) then

         If (a.NE.dble(param(2))) check1=.TRUE.
         If (rmin.NE.dble(param(3))) check1=.TRUE.
         If (rmax.NE.dble(param(4))) check1=.TRUE.
         If (obs.NE.dble(param(5))) check1=.TRUE.
         If (index.NE.dble(param(6))) check1=.TRUE.
         If (rtype.NE.int(param(7))) check1=.TRUE.
         If (mtype.NE.int(param(8))) check1=.TRUE.

      Else if (ftype.EQ.3) then

         If (a.NE.dble(param(2))) check1=.TRUE.
         If (rmax.NE.dble(param(4))) check1=.TRUE.
         If (obs.NE.dble(param(5))) check1=.TRUE.
         If (index.NE.dble(param(6))) check1=.TRUE.
         If (rtype.NE.int(param(7))) check1=.TRUE.
         If (mtype.NE.int(param(8))) check1=.TRUE.

         If (rmin.NE.dble(param(3))) then

            risco = rms(dble(param(2)))

            If ((rmin.NE.dble(param(3))).AND.
     &          (dble(param(3)).GE.risco)) check1=.TRUE.

            If ((rmin.GT.risco).AND.
     &           (dble(param(3)).LT.risco)) check1=.TRUE.

         End If

      End If

      !write(*,*)'Does transfer function require calculation?',check1


      check2=.NOT.(check1)

      If (check1) then


         a = dble(param(2))
         rmin = dble(param(3))
         rmax = dble(param(4))
         obs = dble(param(5))
         index = dble(param(6))
         rtype = int(param(7))
         mtype = int(param(8))
         ftype = int(param(9))

         !write(*,*)'Line parameters',rmin,rmax,obs,index,rtype,mtype


         check4=.FALSE.

         rval = rms(0.998d0)

         If ((ftype.EQ.0).AND.(rtype.EQ.2)
     &        .AND.(rmin.LT.rval)) check4=.TRUE.

         If ((ftype.EQ.1).AND.(rtype.EQ.2)
     &        .AND.(a.GT.0.998d0)) check4=.TRUE.

         If (check4) then

            rtype=0

            write(*,*)'WARNING - Zero Torque Inner Boundary condition
     &is invalid for a = 1.0, resetting radial emissivity type'

         End If


         call lineinit (nmask,n,e,de,nb)

         !write(*,*)'Initialised line'


         check2=.FALSE.

         call obsdet (check2,nmask,n,e,de,r1,r2,rmin,rmax,a,as,obs,
     &               index,rtype,mtype,ftype,nb,ni,ngp,ng,nr,nc,np)

         !write(*,*)'Calculated line for inclination, result',check2


         param(2) = real(a)
         param(3) = real(rmin)
         param(4) = real(rmax)
         param(5) = real(obs)
         param(6) = real(index)
         param(7) = real(rtype)
         param(8) = real(mtype)

         !write(*,*)'Passed back parameters'


      End If


      check4=.FALSE.
      If (check2) check4=.TRUE.
      If (nes.NE.ne) check4=.TRUE.
      If (ec.NE.dble(param(1))) check4=.TRUE.
      If (emin.NE.dble(ear(0))) check4=.TRUE.
      If (emax.NE.dble(ear(ne-1))) check4=.TRUE.

      !write(*,*)'Does line require rebinning?',check4


      If (check4) then


         nes = ne
         emin = dble(ear(0))
         emax = dble(ear(ne-1))
         ec = dble(param(1))

         !write(*,*)'Line parameters',nes,emin,emax,ec


         qfirst=.FALSE.

         call rebinp2 (ne,ec,ear,photar,nmask,n,e,de,nb)

         !write(*,*)'Rebinned line onto XSPEC grid'


      End If


      as = a

      fts=ftype


      If (ccheck) then


         write(buff,'(a,f6.2)')'Spin=',a

         call xwrite(buff,12)

         write(buff,'(a,f6.2)')'Rmin(G)=',r1

         call xwrite(buff,12)

         write(buff,'(a,f6.2)')'Rmax(G)=',r2

         call xwrite(buff,12)


      End If


      return

      end


c *************************************************************************
c *               These are the subroutines for xskdline                  *
c *************************************************************************


      subroutine lineinit (nmask,n,e,de,ne)

      implicit none

      integer ne
      double precision n(ne)
      double precision e(0:ne-1)
      double precision de(0:ne-1)
      logical nmask(0:ne-1)

      integer i


      do i=1,ne-1

         n(i-1)=0.0d0

         e(i-1)=0.0d0

         de(i-1)=0.0d0

         nmask(i-1)=.FALSE.

      end do


      return

      end



      subroutine obsdet (checkout,nmask,n,e,de,r1,r2,rmin,rmax,a,as,obs,
     &                    index,rtype,mtype,ftype,nb,ni,ngp,ng,nr,nc,np)

      implicit none

      integer nb,ni,ngp
      integer ng,nr,nc,np
      integer rtype,mtype,ftype
      double precision r1,r2
      double precision rmin,rmax
      double precision de(0:nb-1)
      double precision a,as,obs,index
      double precision e(0:nb-1),n(nb)
      logical checkout,nmask(0:nb-1)

      integer no,nb1
      parameter (no=2,nb1=100)

      double precision oval(no)
      double precision gb(no,2)
      double precision g(no,0:nb1-1)
      double precision dg(no,0:nb1-1)
      double precision nphot(no,nb1)

      logical omask(no),check1
      logical gmask(no,0:nb1-1)


      call ovaldet (omask,oval,obs,no)

      !write(*,*)'Inclinations to use',oval(1),oval(2)


      check1=.FALSE.

      call olinedet (check1,gmask,nphot,g,dg,omask,oval,r1,r2,rmin,rmax,
     &           a,as,index,rtype,mtype,ftype,no,nb1,ni,ngp,ng,nr,nc,np)

      !write(*,*)'Determined lines for inclinations, result',check1


      checkout=.FALSE.

      If (check1) then


         call olmergep1 (omask,gmask,gb,nphot,g,dg,no,nb1)

         !write(*,*)'Placed lines onto common grid'


         call olmergep2 (checkout,omask,oval,obs,nmask,n,e,de,gmask,
     &                                      gb,nphot,g,dg,nb,no,nb1)

         !write(*,*)'Interpolated between lines, result',checkout


      End If


      return

      end


c *************************************************************************
c *                 These are the subroutines for obsdet                  *
c *************************************************************************



      subroutine ovaldet (omask,oval,obs,no)

      implicit none

      integer no
      double precision oval(no),obs
      logical omask(no)

      double precision o


      o=dble(nint(obs))

      !write(*,*)'Integer form of inclination',o,obs


      If (o.EQ.obs) then

         oval(1)=obs

         oval(2)=0.0d0

         omask(1)=.TRUE.

         omask(2)=.FALSE.

      Else if (o.LT.obs) then

         oval(1)=o

         oval(2)=o+1.0d0

         omask(1)=.TRUE.

         omask(2)=.TRUE.

      Else if (o.GT.obs) then

         oval(2)=o

         oval(1)=o-1.0d0

         omask(1)=.TRUE.

         omask(2)=.TRUE.

      End If


      return

      end



      subroutine olinedet (checkout,nmask,n,e,de,omask,oval,r1,r2,
     &  rmin,rmax,a,as,index,rtype,mtype,ftype,no,nb,ni,ngp,ng,nr,nc,np)

      implicit none

      integer ng,nr,nc,np
      integer no,nb,ni,ngp
      integer rtype,mtype,ftype
      double precision r1,r2
      double precision rmin,rmax
      double precision a,as,index
      logical checkout

      double precision oval(no)
      double precision n(no,nb)
      double precision e(no,0:nb-1)
      double precision de(no,0:nb-1)

      logical omask(no)
      logical nmask(no,0:nb-1)

      integer nb1
      parameter (nb1=100)

      double precision obs
      double precision g(0:nb1-1)
      double precision dg(0:nb1-1)
      double precision nphot(nb1)

      logical gmask(0:nb1-1)
      logical check1,check2

      integer n1,n2
      integer i,j


      n1=0

      n2=0

      do i=1,no


         check1=omask(i)

         If (check1) then


            n1=n1+1

            obs=oval(i)

            check2=.FALSE.

            If (ftype.EQ.0) call lineinit0 (check2,gmask,nphot,g,dg,
     &   r1,r2,rmin,rmax,a,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

            If (ftype.EQ.1) call lineinit1 (check2,gmask,nphot,g,dg,
     &r1,r2,rmin,rmax,a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

            If (ftype.EQ.2) call lineinit2 (check2,gmask,nphot,g,dg,
     &r1,r2,rmin,rmax,a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

            !write(*,*)i,' Calculated line for fit type, result',check1


            If (check2) then


               n2=n2+1

               omask(i)=.TRUE.

               do j=1,nb

                  n(i,j)=nphot(j)

                  e(i,j-1)=g(j-1)

                  de(i,j-1)=dg(j-1)

                  nmask(i,j-1)=gmask(j-1)

               end do

               !write(*,*)i,' Stored line profiles'


            End If


         End If


      end do


      checkout=.FALSE.

      If (n1.EQ.n2) checkout=.TRUE.

      !write(*,*)'Were valid inclinations found?',checkout,n1,n2


      return

      end



      subroutine olmergep1 (omask,nmask,gb,n,e,de,no,nb)

      implicit none

      integer no,nb
      double precision gb(no,2)
      double precision n(no,nb)
      double precision e(no,0:nb-1)
      double precision de(no,0:nb-1)

      logical omask(no)
      logical nmask(no,0:nb-1)

      integer nb1
      parameter (nb1=100)

      double precision gmin,gmax,x,y
      double precision nphot(nb1),dg1
      double precision g(nb1),dg(nb1)

      logical gmask(nb1),check1

      integer i,j


      If (omask(1).AND.omask(2)) then


         gmin=1.0d2

         gmax=-1.0d2

         do i=1,no

            gb(i,1)=e(i,0)

            gb(i,2)=e(i,nb-1)

            gmin=min(gmin,e(i,0))

            gmax=max(gmax,e(i,nb-1))

         end do

         !write(*,*)'Limits on energy grid',gmin,gmax


         dg1=(gmax-gmin)/dble(nb-1)

         !write(*,*)'Bin width',dg1


         do i=1,no


            do j=1,nb

               dg(j)=de(i,j-1)

               gmask(j)=nmask(i,j-1)

               nphot(j)=n(i,j)/de(i,j-1)

               g(j)=e(i,j-1)+0.5d0*de(i,j-1)

            end do

            !write(*,*)i,' Setup interpolation arrays'


            e(i,0)=gmin

            do j=1,nb


               n(i,j)=0.0d0

               de(i,j-1)=dg1

               nmask(i,j-1)=.FALSE.

               If (j.GT.1) e(i,j-1)=e(i,j-2)+dg1

              !write(*,*)i,j,' Initialised global grid'


               y=0.0d0

               check1=.FALSE.

               x=e(i,j-1)+0.5d0*dg1

               call linterpolate (check1,y,x,gmask,g,nphot,nb1)

               !write(*,*)i,j,' Interpolation',x,y,check1

               If (check1) then

                  n(i,j)=y*dg1

                  nmask(i,j-1)=.TRUE.

               End If

               !write(50,*)i,j-1,e(i,j-1),de(i,j-1),n(i,j)



            end do

            !write(50,*)


         end do


      End If


      return

      end




      subroutine olmergep2 (checkout,omask,oval,obs,nmask,n,e,de,gmask,
     &                                          gb,nphot,g,dg,nb,ni,ng)

      implicit none

      integer nb,ni,ng
      double precision gb(ni,2)
      double precision oval(ni),obs
      double precision nphot(ni,ng)
      double precision de(0:nb-1)
      double precision g(ni,0:ng-1)
      double precision dg(ni,0:ng-1)
      double precision e(0:nb-1),n(nb)
      logical omask(ni),gmask(ni,0:nb-1)
      logical checkout,nmask(0:nb-1)

      integer nb1,ni1,ng1
      parameter (nb1=100,ni1=2,ng1=100)

      double precision work1a(nb1)
      double precision work1d(nb1)
      double precision work1e(nb1)
      double precision work1f(nb1)
      double precision gmin,gmax,de1
      double precision work1b(0:nb1-1)
      double precision work1c(0:nb1-1)
      double precision rfrac,risco,rms,x,y
      double precision aval(ni1),rval(ni1)

      logical mask1a(ni1)
      logical mask1b(0:nb1-1)
      logical mask1c(nb1),check1,check2

      integer i,j,nvalid



      rfrac=1.0d0

      If (omask(1).AND.omask(2)) rfrac=(obs-oval(1))/(oval(2)-oval(1))

      !write(*,*)'Fractional distance between inclinations',rfrac


      gmin=g(1,0)

      gmax=g(1,ng-1)

      de1=(gmax-gmin)/dble(nb1-1)

      !write(*,*)'Lower, upper limits of line',gmin,gmax

      !write(*,*)'Bin width',de1


      nvalid=0

      work1b(0)=gmin

      do i=1,nb1


         work1a(i)=0.0d0

         work1c(i-1)=de1

         mask1b(i-1)=.FALSE.

         If (i.GT.1) work1b(i-1)=work1b(i-2)+de1

         If (omask(1).AND.omask(2)) then

            If (gmask(1,i-1).AND.gmask(2,i-1)) then

               nvalid=nvalid+1

               mask1b(i-1)=.TRUE.

               work1a(i)=min(nphot(1,i),nphot(2,i))+
     &           dabs(rfrac*(nphot(2,i)-nphot(1,i)))

            Else

               If (gmask(1,i-1).AND.((work1b(i-1).LT.gb(2,1)).OR.
     &             (work1b(i-1).GT.gb(2,2)))) then

                  nvalid=nvalid+1

                  mask1b(i-1)=.TRUE.

                  work1a(i)=rfrac*nphot(1,i)

               Else if (gmask(2,i-1).AND.((work1b(i-1).LT.gb(1,1)).OR.
     &             (work1b(i-1).GT.gb(1,2)))) then

                  nvalid=nvalid+1

                  mask1b(i-1)=.TRUE.

                  work1a(i)=rfrac*nphot(2,i)

               End If

            End If

            !write(48,*)i,work1b(i-1),nphot(1,i),work1a(i),nphot(2,i),
!     &                 gmask(1,i-1),mask1b(i-1),gmask(2,i-1)

         Else if (omask(1)) then

            If (gmask(1,i-1)) then

               nvalid=nvalid+1

               work1a(i)=nphot(1,i)

               mask1b(i-1)=.TRUE.

            End If

         Else if (omask(2)) then

            If (gmask(2,i-1)) then

               nvalid=nvalid+1

               work1a(i)=nphot(2,i)

               mask1b(i-1)=.TRUE.

            End If

         End If


      end do


      check1=.FALSE.

      If (nvalid.GT.0) check1=.TRUE.

      !write(*,*)'Were valid points found in line?',check1


      nvalid=0

      If (check1) then


         gmin=work1b(0)

         gmax=work1b(nb1-1)

         If (omask(1).AND.omask(2)) then


            gmin=min(gb(1,1),gb(2,1))+dabs(rfrac*(gb(2,1)-gb(1,1)))

            gmax=max(gb(1,1),gb(2,2))+dabs(rfrac*(gb(2,2)-gb(1,2)))


            do i=1,nb1

               mask1c(i)=mask1b(i-1)

               work1f(i)=work1c(i-1)

               work1d(i)=work1a(i)/work1c(i-1)

               work1e(i)=work1b(i-1)+0.5d0*work1c(i-1)

               !write(49,*)work1e(i),work1d(i)*work1c(i-1)

            end do

            !write(49,*)'no no'

            !write(*,*)'Setup interpolation arrays'


         End If

         !write(*,*)'Lower, upper limits on line',gmin,gmax


         de1=(gmax-gmin)/dble(nb1-1)

         !write(*,*)'Bin width',de1


         e(0)=gmin

         do i=1,nb


            n(i)=0.0d0

            de(i-1)=de1

            nmask(i-1)=.FALSE.

            If (i.GT.1) e(i-1)=e(i-2)+de1

            If (omask(1).AND.omask(2)) then


               y=0.0d0

               check2=.FALSE.

               x=e(i-1)+0.5d0*de1

               call linterpolate (check2,y,x,mask1c,work1e,work1d,nb1)

              !write(*,*)i,' Line interpolation',x,y,check2


               If (check2) then

                  n(i)=y*de1

                  nvalid=nvalid+1

                  nmask(i-1)=.TRUE.

               End If


            Else


               If (mask1b(i-1)) then


                  n(i)=work1a(i)

                  nvalid=nvalid+1

                  nmask(i-1)=.TRUE.


               End If


            End If

            !write(49,*)e(i-1),n(i)


         end do

         !write(49,*)'no no'


      End If


      checkout=.FALSE.

      If (nvalid.GT.0) checkout=.TRUE.

      !write(*,*)'Were valid points found in line?',checkout


      return

      end



      subroutine lineinit0 (checkout,nmask,n,e,de,r1,r2,rmin,rmax,
     &        a,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

      implicit none

      integer ngp,ng,nr,nc,np
      integer rtype,mtype,nb,ni
      double precision r1,r2
      double precision rmin,rmax
      double precision de(0:nb-1)
      double precision a,obs,index
      double precision e(0:nb-1),n(nb)
      logical checkout,nmask(0:nb-1)

      integer nb1,ni1,ngp1,ng1,nr1,nc1,np1
      parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5)

      double precision gb(ni1,2)
      double precision g2(ni1,0:nb1-1)
      double precision dg2(ni1,0:nb1-1)
      double precision nphot2(ni1,nb1)
      double precision nphot1(ni1,ngp1,nc1,ng1)
      double precision g1(ni1,ngp1,nc1,0:ng1-1)
      double precision dg1(ni1,ngp1,nc1,0:ng1-1)
      double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1)

      logical gmask2(ni1,0:nb1-1)
      logical gpmask(ni1,ngp1),amask(ni1)
      logical gmask1(ni1,ngp1,nc1,0:ng1-1)
      logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1)

      logical check1,check2

      character*80 buff

      integer ftype


      r1=rmin

      r2=rmax

      call adet (a,rmin)

      !write(*,*)'Spin of black hole',a


      ftype=0

      call tfset (check1,amask,gpmask,tfmask,tf,rmin,rmax,a,obs,
     &                           ftype,ni1,ngp1,ng1,nr1,nc1,np1)

      !write(*,*)'Completed transfer function read, result',check1


      checkout=.FALSE.


      If (check1) then


         call ldet (check2,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask,tf,
     &                        index,a,rtype,mtype,ni,ngp,ng,nr,nc,np)

         !write(*,*)'Determined line profiles, result',check2


         If (check2) then


            call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask,
     &                gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1)

            !write(*,*)'Merged line profiles onto common energy grid'


            call lmergep2 (checkout,nmask,n,e,de,amask,
     &            gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1)

            !write(*,*)'Determined line profile for spin'


         End If


      End If


      return

      end



      subroutine lineinit1 (checkout,nmask,n,e,de,r1,r2,rmin,rmax,
     &     a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

      implicit none

      integer ngp,ng,nr,nc,np
      integer rtype,mtype,nb,ni
      double precision r1,r2
      double precision rmin,rmax
      double precision de(0:nb-1)
      double precision a,as,obs,index
      double precision e(0:nb-1),n(nb)
      logical checkout,nmask(0:nb-1)

      integer nb1,ni1,ngp1,ng1,nr1,nc1,np1
      parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5)

      double precision gb(ni1,2)
      double precision risco,rms
      double precision g2(ni1,0:nb1-1)
      double precision dg2(ni1,0:nb1-1)
      double precision nphot2(ni1,nb1)
      double precision nphot1(ni1,ngp1,nc1,ng1)
      double precision g1(ni1,ngp1,nc1,0:ng1-1)
      double precision dg1(ni1,ngp1,nc1,0:ng1-1)
      double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1)

      logical gmask2(ni1,0:nb1-1)
      logical gpmask(ni1,ngp1),amask(ni1)
      logical gmask1(ni1,ngp1,nc1,0:ng1-1)
      logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1)

      logical check1,check2,check4

      character*80 buff

      integer ftype


      check1=.FALSE.

      call rconvert (check1,r1,r2,rmin,rmax,a)

      !write(*,*)'Converted radii to gravitational units',r1,r2,check1


      If (check1) then


         risco = rms(a)

         If (r1.LT.risco) r1 = risco

         !write(*,*)'New r_{ms}',risco

         !write(*,*)'Revised inner edge of disc',r1

         !write(*,*)'Spin of black hole',a

         !write(*,*)'Inner edge of disc, Rin(G)',r1

         !write(*,*)'Outer edge of disc, Rout(G)',r2


         ftype=1

         call tfset (check2,amask,gpmask,tfmask,tf,r1,rmax,a,obs,
     &                            ftype,ni1,ngp1,ng1,nr1,nc1,np1)

         !write(*,*)'Completed transfer function read, result',check2


         checkout=.FALSE.

         If (check2) then


            call ldet (check4,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask,
     &                       tf,index,a,rtype,mtype,ni,ngp,ng,nr,nc,np)

            !write(*,*)'Determined line profiles, result',check4


            If (check4) then


               call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask,
     &                   gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1)

               !write(*,*)'Merged line profiles onto common energy grid'


               call lmergep2 (checkout,nmask,n,e,de,amask,
     &               gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1)

               !write(*,*)'Determined line profile for spin'


            End If


         End If


      End If


      return

      end


      subroutine lineinit2 (checkout,nmask,n,e,de,r1,r2,rmin,rmax,
     &     a,as,obs,index,rtype,mtype,nb,ni,ngp,ng,nr,nc,np)

      implicit none

      integer ngp,ng,nr,nc,np
      integer rtype,mtype,nb,ni
      double precision r1,r2
      double precision rmin,rmax
      double precision de(0:nb-1)
      double precision a,as,obs,index
      double precision e(0:nb-1),n(nb)
      logical checkout,nmask(0:nb-1)

      integer nb1,ni1,ngp1,ng1,nr1,nc1,np1
      parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5)

      double precision gb(ni1,2)
      double precision risco,rms
      double precision g2(ni1,0:nb1-1)
      double precision dg2(ni1,0:nb1-1)
      double precision nphot2(ni1,nb1)
      double precision aval(ni1),rval(ni1)
      double precision nphot1(ni1,ngp1,nc1,ng1)
      double precision g1(ni1,ngp1,nc1,0:ng1-1)
      double precision dg1(ni1,ngp1,nc1,0:ng1-1)
      double precision tf(ni1,ngp1,ng1+1,nr1+1,nc1,np1)

      logical gmask2(ni1,0:nb1-1)
      logical gpmask(ni1,ngp1),amask(ni1)
      logical gmask1(ni1,ngp1,nc1,0:ng1-1)
      logical tfmask(ni1,ngp1,ng1+1,nr1+1,nc1)

      logical check1,check2,check4

      integer ftype


      call aclose (amask,rval,aval,a,ni1,ni1)

      write(*,*)'Closest values in a',aval(1),aval(2)


      check1=amask(1).AND.amask(2)

      If (check1) then

         r1=dabs(aval(1)-a)

         r2=dabs(aval(2)-a)

         If (r1.LE.r2) amask(2)=.FALSE.

         If (r1.GT.r2) amask(1)=.FALSE.

         If (amask(1)) a=aval(1)

         If (amask(2)) a=aval(2)

         If (amask(1)) as=aval(1)

         If (amask(2)) as=aval(2)

         write(*,*)'Spin parameters',a,aval(1),aval(2),amask(1),amask(2)

      End If


      check1=.FALSE.

      call rconvert (check1,r1,r2,rmin,rmax,a)

      write(*,*)'Converted radii to gravitational units',r1,r2,check1


      If (check1) then


         risco = rms(a)

         If (r1.LT.risco) r1 = risco

         write(*,*)'New r_{ms}',risco

         write(*,*)'Revised inner edge of disc',r1


         ftype=2

         call tfset (check2,amask,gpmask,tfmask,tf,r1,rmax,a,obs,
     &                            ftype,ni1,ngp1,ng1,nr1,nc1,np1)

         write(*,*)'Completed transfer function read, result',check2


         checkout=.FALSE.

         If (check2) then


            call ldet (check4,gmask1,nphot1,g1,dg1,amask,gpmask,tfmask,
     &                       tf,index,a,rtype,mtype,ni,ngp,ng,nr,nc,np)

            write(*,*)'Determined line profiles, result',check4


            If (check4) then


               call lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask,
     &                   gmask1,nphot1,g1,dg1,nb1,ni1,ngp1,nc1,ng1)

               write(*,*)'Merged line profiles onto common energy grid'


               call lmergep2 (checkout,nmask,n,e,de,amask,
     &               gmask2,gb,nphot2,g2,dg2,a,nb,ni1,nb1)

               write(*,*)'Determined line profile for spin'


            End If


         End If


      End If


      return

      end



c *************************************************************************
c *               These are the subroutines for lineinit                  *
c *************************************************************************



      subroutine adet (a,rmin)

      implicit none

      double precision a,rmin

      double precision a1,a2
      double precision aval(101)
      double precision rval(101)
      double precision scale,rms

      logical check1,mask1(101)

      integer i,n


      n=101

      a1=1.0d0

      a2=0.0d0

      scale=1.0d-2

      call gridinit (aval,a1,a2,scale,n)

      !write(*,*)'Setup spin grid'


      do i=1,n

         mask1(i)=.TRUE.

         rval(i)=rms(aval(i))

      end do

      !write(*,*)'Initialised r_{ms} grid'


      check1=.FALSE.

      call linterpolate (check1,a,rmin,mask1,rval,aval,n)

      !write(*,*)'Interpolated value of a',rmin,a,check1


      return

      end



      subroutine rconvert (checkout,r1,r2,rmin,rmax,a)

      implicit none

      double precision r1,r2,rmin,rmax,a
      logical checkout

      integer n
      parameter (n=1000)

      double precision xmin,xmax
      double precision risco,rms,rout
      double precision x(n),y(n),scale

      logical mask1(n),check1,check2

      integer i


      xmin=1.0d0

      xmax=1.0d2

      rout=1.0d2

      scale=1.0d-3

      risco = rms(a)

      call gridinit (x,xmin,xmax,scale,n)

      call gridinit (y,risco,rout,scale,n)

      !write(*,*)'Initialised radial conversion grids'


      do i=1,n

         mask1(i)=.TRUE.

         !write(60,*)i,x(i),y(i),mask1(i)

      end do

      !write(*,*)'Initialised interpolation mask'


      check1=.FALSE.

      check2=.FALSE.

      call linterpolate (check1,r1,rmin,mask1,x,y,n)

      call linterpolate (check2,r2,rmax,mask1,x,y,n)

      !write(*,*)'Gravitational radii',r1,r2,check1,check2


      checkout=check1.AND.check2

      !write(*,*)'Were conversions succesfull?',checkout


      return

      end



      subroutine tfset (checkout,amask,gpmask,tfmask,tf,rmin,rmax,a,obs,
     &                                         ftype,ni,ngp,ng,nr,nc,np)

      implicit none

      integer ftype
      integer ni,ngp,ng,nr,nc,np
      double precision rmin,rmax,a,obs
      double precision tf(ni,ngp,ng+1,nr+1,nc,np)
      logical tfmask(ni,ngp,ng+1,nr+1,nc)
      logical gpmask(ni,ngp),amask(ni)
      logical checkout

      integer ni1,ngp1,ng1,nr1,nc1,np1
      parameter (ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5)

      double precision aval(ni1),rval(ni1),r1,r2,a1
      double precision work5(ngp1,ng1+1,nr1+1,nc1,np1)

      logical check1,check2,mask1(ngp1)
      logical mask4(ngp1,ng1+1,nr1+1,nc1)

      integer n1,n2
      integer i,j,k,l,m,n


      call aclose (amask,rval,aval,a,ni,ni1)

      !write(*,*)'Closest values in a',aval(1),aval(2)


      check1=amask(1).AND.amask(2)

      If ((ftype.EQ.2).AND.check1) then

         r1=dabs(aval(1)-a)

         r2=dabs(aval(2)-a)

         If (r1.LE.r2) amask(2)=.FALSE.

         If (r1.GT.r2) amask(1)=.FALSE.

         !write(*,*)'Spin parameters',aval(1),aval(2),amask(1),amask(2)

      End If


      n1=0

      n2=0

      do i=1,ni


         check1=amask(i)

         If (check1) then


            r1=rmin

            r2=rmax

            n1=n1+1

            a1=aval(i)

            check2=.FALSE.

            If (ftype.EQ.0) r1=rval(i)

            call tfread (check2,mask1,mask4,work5,r1,r2,a1,obs,
     &                                    ngp1,nc1,ng1,nr1,np1)

            !write(*,*)i,' Read in transfer function, result',check2


            amask(i)=check2

            If (amask(i)) then


               call rminset (mask1,mask4,work5,r1,a1,
     &                          ngp1,nc1,ng1,nr1,np1)

               !write(*,*)i,' Regridded inner edge'


               call rmaxset (mask1,mask4,work5,r2,a1,
     &                          ngp1,nc1,ng1,nr1,np1)

               !write(*,*)i,' Regridded outer edge'


               n2=n2+1

               do j=1,ngp

                  gpmask(i,j)=mask1(j)

                  If (mask1(j)) then

                     do k=1,ng+1

                        do l=1,nr+1

                           do m=1,nc

                              tfmask(i,j,k,l,m)=mask4(j,k,l,m)

                              If (mask4(j,k,l,m)) then

                                 do n=1,np

                                    tf(i,j,k,l,m,n)=work5(j,k,l,m,n)

                                 end do

                              End If

                           end do

                        end do

                     end do

                  End If

               end do

               !write(*,*)i,' Placed transfer function in storage arrays'


            End If


         End If

      end do


      checkout=.FALSE.

      If (n1.EQ.n2) checkout=.TRUE.

      !write(*,*)'Were transfer functions read in successfully?',checkout


      return

      end



      subroutine ldet (checkout,gmask,nphot,g,dg,amask,gpmask,tfmask,tf,
     &                           index,a,rtype,mtype,ni,ngp,ng,nr,nc,np)

      implicit none

      integer rtype,mtype
      integer ni,ngp,ng,nr,nc,np
      double precision rmin,rmax
      double precision a,obs,index
      double precision nphot(ni,ngp,nc,ng)
      double precision g(ni,ngp,nc,0:ng-1)
      double precision dg(ni,ngp,nc,0:ng-1)
      double precision tf(ni,ngp,ng+1,nr+1,nc,np)
      logical tfmask(ni,ngp,ng+1,nr+1,nc)
      logical gpmask(ni,ngp),amask(ni)
      logical gmask(ni,ngp,nc,0:ng-1)
      logical checkout

      integer ni1,ngp1,ng1,nr1,nc1,np1
      parameter (ni1=2,ngp1=5,nc1=2,ng1=100,nr1=100,np1=5)

      double precision work3a(ngp1,nc1,ng1)
      double precision work3b(ngp1,nc1,0:ng1-1)
      double precision work3c(ngp1,nc1,0:ng1-1)
      double precision work5(ngp1,ng1+1,nr1+1,nc1,np1)

      logical mask1(ngp1),check1
      logical mask3(ngp1,nc1,0:ng1-1)
      logical mask4(ngp1,ng1+1,nr1+1,nc1)

      integer n1,n2
      integer i,j,k,l,m,n


      n1=0

      n2=0

      do i=1,ni


         check1=amask(i)

         If (check1) then


            n1=n1+1

            do j=1,ngp

               mask1(j)=gpmask(i,j)

               do k=1,ng+1

                  do l=1,nr+1

                     do m=1,nc

                        mask4(j,k,l,m)=tfmask(i,j,k,l,m)

                        do n=1,np

                           work5(j,k,l,m,n)=tf(i,j,k,l,m,n)

                        end do

                     end do

                  end do

               end do

            end do

            !write(*,*)i,' Extracted transfer function'


            amask(i)=.FALSE.

            call lcalc (amask(i),work3b,work3a,work3c,
     &                  mask3,mask1,mask4,work5,
     &                  index,a,rtype,mtype,ngp1,nc1,ng1,nr1,np1)

            !write(*,*)i,' Calculated line profiles, result',amask(i)


            If (amask(i)) then


               n2=n2+1

               do j=1,ngp

                  If (gpmask(i,j)) then

                     do k=1,nc

                        do l=1,ng

                           nphot(i,j,k,l)=work3a(j,k,l)
                           g(i,j,k,l-1)=work3b(j,k,l-1)
                           dg(i,j,k,l-1)=work3c(j,k,l-1)
                           gmask(i,j,k,l-1)=mask3(j,k,l-1)
                           !write(44,*)i,j,k,l,dg(i,j,k,l-1)

                        end do

                      end do

                   End If

                end do

                !write(*,*)i,' Stored line profile'



            End If


         End If


      end do


      checkout=.FALSE.

      If (n1.EQ.n2) checkout=.TRUE.

      !write(*,*)'Were line profiles calculated successfully?',checkout


      return

      end



      subroutine lmergep1 (gmask2,gb,nphot2,g2,dg2,amask,gpmask,
     &                     gmask1,nphot1,g1,dg1,nb,ni,ngp,nc,ng)

      implicit none

      integer nb,ni,ngp,nc,ng
      double precision gb(ni,2)
      double precision g2(ni,0:nb-1)
      double precision dg2(ni,0:nb-1)
      double precision nphot2(ni,nb)
      double precision nphot1(ni,ngp,nc,ng)
      double precision g1(ni,ngp,nc,0:ng-1)
      double precision dg1(ni,ngp,nc,0:ng-1)
      logical gmask1(ni,ngp,nc,0:ng-1)
      logical gpmask(ni,ngp),amask(ni)
      logical gmask2(ni,0:nb-1)

      integer nb1,ni1,ngp1,nc1,ng1
      parameter (nb1=100,ni1=2,ngp1=5,nc1=2,ng1=100)

      double precision gmin,gmax
      double precision work1a(nb1)
      double precision work1b(0:nb-1)
      double precision work1c(0:nb-1)
      double precision work3a(ngp1,nc1,ng1)
      double precision work3b(ngp1,nc1,0:ng1-1)
      double precision work3c(ngp1,nc1,0:ng1-1)

      logical mask1a(ngp)
      logical mask1b(0:nb-1)
      logical mask3(ngp1,nc1,0:ng1-1)

      integer i,j,k,l


      gmin=1.0d2

      gmax=-1.0d2

      do i=1,ni

         gb(i,1)=0.0d0

         gb(i,2)=0.0d0

         If (amask(i)) then

            gb(i,1)=1.0d2

            gb(i,2)=-1.0d2

            do j=1,ngp

               If (gpmask(i,j)) then

                  do k=1,nc

                     gmin=min(gmin,g1(i,j,k,0))

                     gmax=max(gmax,g1(i,j,k,ng-1))

                     gb(i,1)=min(gb(i,1),g1(i,j,k,0))

                     gb(i,2)=max(gb(i,2),g1(i,j,k,ng-1))

                  end do

               End If

            end do

         End If

      end do

      !write(*,*)'Lower, upper limits of line',gmin,gmax


      do i=1,ni


         If (amask(i)) then


            do j=1,ngp

               mask1a(j)=gpmask(i,j)

               If (gpmask(i,j)) then

                  do k=1,nc

                     do l=1,ng

                        work3a(j,k,l)=nphot1(i,j,k,l)
                        work3b(j,k,l-1)=g1(i,j,k,l-1)
                        work3c(j,k,l-1)=dg1(i,j,k,l-1)
                        mask3(j,k,l-1)=gmask1(i,j,k,l-1)

                        !write(46,*)i,j,k,l,work3b(j,k,l-1),
!     &                  work3a(j,k,l),work3c(j,k,l-1),
!     &                  mask3(j,k,l-1),mask1a(j)

                     end do

                  end do

               End If

            end do

            !write(*,*)i,' Extracted line profile'


            call rebinp1 (mask1b,work1a,work1b,work1c,gmin,gmax,
     &                        mask1a,mask3,work3a,work3c,work3b,
     &                                         nb1,ngp1,ng1,nc1)

            !write(*,*)i,' Rebinned line profiles'


            call lsmooth (mask1b,work1a,work1b,work1c,mask1a,
     &                         mask3,work3b,nb1,ngp1,ng1,nc1)

            !write(*,*)i,' Smoothed line boundaries'


            do j=1,nb

               nphot2(i,j)=work1a(j)

               g2(i,j-1)=work1b(j-1)

               dg2(i,j-1)=work1c(j-1)

               gmask2(i,j-1)=mask1b(j-1)

               !write(40,*)i,j,work1b(j-1),work1c(j-1),
!     &                    work1a(j),mask1b(j-1)

               !write(41,*)i,j,g2(i,j-1),dg2(i,j-1),
!     &                    nphot2(i,j),gmask2(i,j-1)

!               If (gmask2(i,j-1)) !write(42,*)g2(i,j-1),nphot2(i,j)

            end do

            !write(*,*)i,' Stored line profiles'

            !write(40,*)

            !write(41,*)

            !write(42,*)'no no'


         End If

      end do


      return

      end



      subroutine lmergep2 (checkout,nmask,n,e,de,amask,
     &                  gmask,gb,nphot,g,dg,a,nb,ni,ng)

      implicit none

      integer nb,ni,ng
      double precision gb(ni,2)
      double precision nphot(ni,ng)
      double precision de(0:nb-1),a
      double precision g(ni,0:ng-1)
      double precision dg(ni,0:ng-1)
      double precision e(0:nb-1),n(nb)
      logical amask(ni),gmask(ni,0:nb-1)
      logical checkout,nmask(0:nb-1)

      integer nb1,ni1,ng1
      parameter (nb1=100,ni1=2,ng1=100)

      double precision work1a(nb1)
      double precision work1d(nb1)
      double precision work1e(nb1)
      double precision work1f(nb1)
      double precision gmin,gmax,de1
      double precision work1b(0:nb1-1)
      double precision work1c(0:nb1-1)
      double precision rfrac,risco,rms,x,y
      double precision aval(ni1),rval(ni1)

      logical mask1a(ni1)
      logical mask1b(0:nb1-1)
      logical mask1c(nb1),check1,check2

      integer i,nvalid


      risco=rms(a)

      !write(*,*)'r_{ms}(a)',risco


      call aclose (mask1a,rval,aval,a,ni1,ni1)

      !write(*,*)'Closest values of r_{ms}(a)',rval(1),rval(2)


      rfrac=1.0d0

      If (amask(1).AND.amask(2)) rfrac=(risco-rval(1))/(rval(2)-rval(1))

      !write(*,*)'Fractional distance between points in r_{ms}',rfrac


      gmin=g(1,0)

      gmax=g(1,ng-1)

      de1=(gmax-gmin)/dble(nb1-1)

      !write(*,*)'Lower, upper limits of line',gmin,gmax

      !write(*,*)'Bin width',de1


      nvalid=0

      work1b(0)=gmin

      do i=1,nb1


         work1a(i)=0.0d0

         work1c(i-1)=de1

         mask1b(i-1)=.FALSE.

         If (i.GT.1) work1b(i-1)=work1b(i-2)+de1

         If (amask(1).AND.amask(2)) then

            If (gmask(1,i-1).AND.gmask(2,i-1)) then

               nvalid=nvalid+1

               mask1b(i-1)=.TRUE.

               work1a(i)=min(nphot(1,i),nphot(2,i))+
     &           dabs(rfrac*(nphot(2,i)-nphot(1,i)))

            Else

               If (gmask(1,i-1).AND.((work1b(i-1).LT.gb(2,1)).OR.
     &             (work1b(i-1).GT.gb(2,2)))) then

                  nvalid=nvalid+1

                  mask1b(i-1)=.TRUE.

                  work1a(i)=rfrac*nphot(1,i)

               Else if (gmask(2,i-1).AND.((work1b(i-1).LT.gb(1,1)).OR.
     &             (work1b(i-1).GT.gb(1,2)))) then

                  nvalid=nvalid+1

                  mask1b(i-1)=.TRUE.

                  work1a(i)=rfrac*nphot(2,i)

               End If

            End If

            !write(48,*)i,work1b(i-1),nphot(1,i),work1a(i),nphot(2,i),
!     &                 gmask(1,i-1),mask1b(i-1),gmask(2,i-1)

         Else if (amask(1)) then

            If (gmask(1,i-1)) then

               nvalid=nvalid+1

               work1a(i)=nphot(1,i)

               mask1b(i-1)=.TRUE.

            End If

         Else if (amask(2)) then

            If (gmask(2,i-1)) then

               nvalid=nvalid+1

               work1a(i)=nphot(2,i)

               mask1b(i-1)=.TRUE.

            End If

         End If


      end do


      check1=.FALSE.

      If (nvalid.GT.0) check1=.TRUE.

      !write(*,*)'Were valid points found in line?',check1


      nvalid=0

      If (check1) then


         gmin=work1b(0)

         gmax=work1b(nb1-1)

         If (amask(1).AND.amask(2)) then


            gmin=min(gb(1,1),gb(2,1))+dabs(rfrac*(gb(2,1)-gb(1,1)))

            gmax=max(gb(1,1),gb(2,2))+dabs(rfrac*(gb(2,2)-gb(1,2)))


            do i=1,nb1

               mask1c(i)=mask1b(i-1)

               work1f(i)=work1c(i-1)

               work1d(i)=work1a(i)/work1c(i-1)

               work1e(i)=work1b(i-1)+0.5d0*work1c(i-1)

            end do

            !write(*,*)'Setup interpolation arrays'


         End If

         !write(*,*)'Lower, upper limits on line',gmin,gmax


         de1=(gmax-gmin)/dble(nb1-1)

         !write(*,*)'Bin width',de1


         e(0)=gmin

         do i=1,nb


            n(i)=0.0d0

            de(i-1)=de1

            nmask(i-1)=.FALSE.

            If (i.GT.1) e(i-1)=e(i-2)+de1

            If (amask(1).AND.amask(2)) then


               y=0.0d0

               check2=.FALSE.

               x=e(i-1)+0.5d0*de1

               call linterpolate (check2,y,x,mask1c,work1e,work1d,nb1)

               !write(*,*)i,' Line interpolation',x,y,check2


               If (check2) then

                  n(i)=y*de1

                  nvalid=nvalid+1

                  nmask(i-1)=.TRUE.

               End If


            Else


               If (mask1b(i-1)) then


                  n(i)=work1a(i)

                  nvalid=nvalid+1

                  nmask(i-1)=.TRUE.


               End If


            End If

            !write(*,*)i,e(i-1),de(i-1),n(i),nmask(i-1)


         end do


      End If


      checkout=.FALSE.

      If (nvalid.GT.0) checkout=.TRUE.

      !write(*,*)'Were valid points found in line?',checkout


      return

      end



      subroutine aclose (amask,rval,aval,a,ni,ni1)

      implicit none

      integer ni,ni1
      double precision rval(ni1),aval(ni1),a
      logical amask(ni)

      double precision aquant(12),rms,tol

      integer i


      do i=1,ni

         aval(i)=0.0d0

         rval(i)=0.0d0

         amask(i)=.FALSE.

      end do

      !write(*,*)'Initialised closest spin arrays'


      do i=1,12

         If (i.EQ.1) aquant(i)=0.0d0

         If (i.EQ.2) aquant(i)=0.1501819899670580d0

         If (i.EQ.3) aquant(i)=0.2940047205828976d0

         If (i.EQ.4) aquant(i)=0.4305113630911228d0

         If (i.EQ.5) aquant(i)=0.5584815599030040d0

         If (i.EQ.6) aquant(i)=0.6763195721910966d0

         If (i.EQ.7) aquant(i)=0.7818758452482436d0

         If (i.EQ.8) aquant(i)=0.8721520262014855d0

         If (i.EQ.9) aquant(i)=0.9428090365807393d0

         If (i.EQ.10) aquant(i)=0.9874884803808748d0

         If (i.EQ.11) aquant(i)=0.998d0

         If (i.EQ.12) aquant(i)=1.0d0

      end do

      !write(*,*)'Initialised transfer function spins'


      tol=1.0d-5

      do i=1,12

         If (dabs(a-aquant(i)).LT.tol) then

            amask(1)=.TRUE.

            aval(1)=aquant(i)

            rval(1)=rms(aval(1))

            !write(*,*)'Closest spin',aval(1)

            exit

         Else if (i.LT.12) then

            If ((aquant(i).LT.a).AND.(aquant(i+1).GT.a)) then


               If (dabs(a-aquant(i+1)).LT.tol) then

                  amask(1)=.TRUE.

                  aval(1)=aquant(i+1)

                  rval(1)=rms(aval(1))

                  !write(*,*)'Closest spin',aval(1)

               Else

                  amask(1)=.TRUE.

                  amask(2)=.TRUE.

                  aval(1)=aquant(i)

                  aval(2)=aquant(i+1)

                  rval(1)=rms(aval(1))

                  rval(2)=rms(aval(2))

                  !write(*,*)'Closest spins',aval(1),aval(2)

               End If


               exit

            End If

         End If

      end do


      return

      end



c *************************************************************************
c *                  These are the generic subroutines                    *
c *************************************************************************



      subroutine tfread (checkout,gpmask,tfmask,tf,
     &                   rmin,rmax,a,obs,ngp,nc,ng,nr,np)

      implicit none

      integer ngp,nc,ng,nr,np
      logical checkout,gpmask(ngp)
      logical tfmask(ngp,ng+1,nr+1,nc)
      double precision rmin,rmax,a,obs
      double precision tf(ngp,ng+1,nr+1,nc,np)

      integer ng1,nr1,nc1
      parameter (nc1=2,ng1=100,nr1=100)

      double precision rms,r1,r2,tol
      double precision v(ng1+1,nr1+1,nc1)
      double precision w(ng1+1,nr1+1,nc1)
      double precision x(ng1+1,nr1+1,nc1)
      double precision y(ng1+1,nr1+1,nc1)
      double precision z(ng1+1,nr1+1,nc1)

      character*256 path

      logical check1
      logical check2
      logical check4
      logical m3(ng1+1,nr1+1,nc1)

      integer nvalid,loc
      integer start,status
      integer i,j,k,l


      start=1

      nvalid = 0

      tol = 1.0d-6

      gpmask(1)=.FALSE.

      If (dabs(a).LT.tol) then

         start=2

         gpmask(1)=.FALSE.

      End If

      !write(*,*)'Line paramaters',rmin,rmax,a,obs

      do i=start,ngp


         r1 = rms(a)

         If (i.EQ.1) r2 = 6.0d0

         If (i.EQ.2) r1 = 6.0d0

         If (i.EQ.2) r2 = 1.0d1

         If (i.EQ.3) r1 = 1.0d1

         If (i.EQ.3) r2 = 2.0d1

         If (i.EQ.4) r1 = 2.0d1

         If (i.EQ.4) r2 = 5.0d1

         If (i.EQ.5) r1 = 5.0d1

         If (i.EQ.5) r2 = 1.0d2

         !write(*,*)i,' Lower limits of radial grouping',r1,rmin,r2

         !write(*,*)i,' Upper limits of radial grouping',r1,rmax,r2


         check1=.FALSE.

         check2=.FALSE.

         check4=.FALSE.

         If ((rmin.LT.r1).AND.(rmin.LT.r2).AND.
     &       (rmax.GT.r1).AND.(rmax.GT.r2)) check1=.TRUE.

         If ((rmin.GE.r1).AND.(rmin.LT.r2)) check2=.TRUE.

         If ((rmax.GT.r1).AND.(rmax.LE.r2)) check4=.TRUE.

         If (check2.OR.check4) check1=.TRUE.

         !write(*,*)i,' Does group contribute?',check1,check2,check4


         gpmask(i)=.FALSE.

         If (check1) then


            call pathinit (path,r1,r2,a,obs)

c            write(*,*)i,' Path to file ',path


            check2=.FALSE.

            call tfopen (check2,m3,v,w,x,y,z,path,ng1,nr1,nc1)

c            write(*,*)i,' Read in transfer function'


            If (check2) then


               gpmask(i)=.TRUE.

               do j=1,ng+1

                  do k=1,nr+1

                     do l=1,nc

                        tf(i,j,k,l,1)=v(j,k,l)

                        tf(i,j,k,l,2)=w(j,k,l)

                        tf(i,j,k,l,3)=x(j,k,l)

                        tf(i,j,k,l,4)=y(j,k,l)

                        tf(i,j,k,l,5)=z(j,k,l)

                        tfmask(i,j,k,l)=m3(j,k,l)

                     end do

                  end do

               end do


            End If

            !write(*,*)i,' Read in transfer function, result',gpmask(i)


         End If


         If (gpmask(i)) nvalid = nvalid + 1


      end do


      checkout=.FALSE.

      If (nvalid.GT.0) checkout=.TRUE.

      !write(*,*)'Was transfer function read succesfull?',checkout


      return

      end



      subroutine rminset (gpmask,tfmask,tf,rmin,a,ngp,nc,ng,nr,np)

      implicit none

      integer ngp,nc,ng,nr,np
      double precision rmin,a
      logical checkout,gpmask(ngp)
      logical tfmask(ngp,ng+1,nr+1,nc)
      double precision tf(ngp,ng+1,nr+1,nc,np)

      integer nr1
      parameter (nr1=101)

      double precision r1,r2,scale,rms
      double precision r(nr1),m(nr1)
      double precision x(nr1),y(nr1)
      double precision z(nr1),w(3)

      logical mask1b(3)
      logical mask1a(nr1)
      logical check1,check2
      logical check4,check5
      logical check6,check7

      integer i,j,k,l


      do i=1,ngp


         If (gpmask(i)) then


            r1 = rms(a)

            If (i.EQ.1) r2 = 6.0d0

            If (i.EQ.2) r1 = 6.0d0

            If (i.EQ.2) r2 = 1.0d1

            If (i.EQ.3) r1 = 1.0d1

            If (i.EQ.3) r2 = 2.0d1

            If (i.EQ.4) r1 = 2.0d1

            If (i.EQ.4) r2 = 5.0d1

            If (i.EQ.5) r1 = 5.0d1

            If (i.EQ.5) r2 = 1.0d2

            !write(*,*)i,' Lower,upper limits of radial grouping',r1,r2


            check1=.FALSE.

            If (rmin.GT.r1) check1=.TRUE.

            !write(*,*)i,' Does grid need reshaping?',r1,rmin,r2,check1


            If (check1) then


               do j=1,nc

                  do k=1,ng


                     r1=0.0d0

                     check2=.FALSE.

                     do l=1,nr+1

                        If (tfmask(i,k,l,j)) then

                           check2=.TRUE.

                           r1=tf(i,k,l,j,2)

                           exit

                        End If

                     end do

                     !write(*,*)i,j,k,' Inner boundary',r1,check2


                     r2=0.0d0

                     check4=.FALSE.

                     do l=nr+1,1,-1

                        If (tfmask(i,k,l,j)) then

                           check4=.TRUE.

                           r2=tf(i,k,l,j,2)

                           exit

                        End If

                     end do

                     !write(*,*)i,j,k,' Outer boundary',r2,check4

                     !write(50,*)i,j,k,r1,r2,check2,check4


                     check5=.FALSE.

                     check6=.FALSE.

                     check7=.FALSE.

                     If (check2.AND.check4) then

                        If ((r1.GT.rmin).AND.(r2.GT.rmin)) check5=.TRUE.

                        If ((r1.LT.rmin).AND.(r2.GT.rmin)) check6=.TRUE.

                        If ((r1.LT.rmin).AND.(r2.LT.rmin)) check7=.TRUE.

                     End If

                     !write(*,*)i,j,k,' Results of boundary tests',
!     &                                       check5,check6,check7

                     !write(51,*)i,j,k,r1,rmin,r2,check5,check6,check7


                     If (check6) then


                        do l=1,nr1

                           r(l)=tf(i,k,l,j,2)

                           m(l)=tf(i,k,l,j,3)

                           x(l)=tf(i,k,l,j,4)

                           y(l)=tf(i,k,l,j,5)

                           mask1a(l)=tfmask(i,k,l,j)

                        end do

                        !write(*,*)i,j,k,' Setup interpolation arrays'


                        scale=1.0d-2

                        call gridinit (z,rmin,r2,scale,nr1)

                        !write(*,*)i,j,k,'Initialised new radial grid'


                        do l=1,nr1

                           w(1)=0.0d0

                           w(2)=0.0d0

                           w(3)=0.0d0

                           mask1b(1)=.FALSE.

                           mask1b(2)=.FALSE.

                           mask1b(3)=.FALSE.

                           tf(i,k,l,j,2)=0.0d0

                           tf(i,k,l,j,3)=0.0d0

                           tf(i,k,l,j,4)=0.0d0

                           tf(i,k,l,j,5)=0.0d0

                           tfmask(i,k,l,j)=.FALSE.

                           call linterpolate (mask1b(1),w(1),z(l),
     &                                        mask1a,r,m,nr1)

                           call linterpolate (mask1b(2),w(2),z(l),
     &                                        mask1a,r,x,nr1)

                           call linterpolate (mask1b(3),w(3),z(l),
     &                                        mask1a,r,y,nr1)

                           If (mask1b(1).AND.mask1b(2)
     &                         .AND.mask1b(3)) then

                              tfmask(i,k,l,j)=.TRUE.

                              tf(i,k,l,j,2)=z(l)

                              tf(i,k,l,j,3)=w(1)

                              tf(i,k,l,j,4)=w(2)

                              tf(i,k,l,j,5)=w(3)

                           End If


                        end do


                     Else if (check7) then


                        do l=1,nr+1

                           tf(i,k,l,j,2)=0.0d0

                           tf(i,k,l,j,3)=0.0d0

                           tf(i,k,l,j,4)=0.0d0

                           tf(i,k,l,j,5)=0.0d0

                           tfmask(i,k,l,j)=.FALSE.

                        end do


                     End If


                  end do

               end do


            End If


            exit


         End If


      end do


      return

      end



      subroutine rmaxset (gpmask,tfmask,tf,rmax,a,ngp,nc,ng,nr,np)

      implicit none

      integer ngp,nc,ng,nr,np
      double precision rmax,a
      logical checkout,gpmask(ngp)
      logical tfmask(ngp,ng+1,nr+1,nc)
      double precision tf(ngp,ng+1,nr+1,nc,np)

      integer nr1
      parameter (nr1=101)

      double precision r1,r2,scale,rms
      double precision r(nr1),m(nr1)
      double precision x(nr1),y(nr1)
      double precision z(nr1),w(3)

      logical mask1b(3)
      logical mask1a(nr1)
      logical check1,check2
      logical check4,check5
      logical check6,check7

      integer i,j,k,l


      do i=ngp,1,-1


         If (gpmask(i)) then


            r1 = rms(a)

            If (i.EQ.1) r2 = 6.0d0

            If (i.EQ.2) r1 = 6.0d0

            If (i.EQ.2) r2 = 1.0d1

            If (i.EQ.3) r1 = 1.0d1

            If (i.EQ.3) r2 = 2.0d1

            If (i.EQ.4) r1 = 2.0d1

            If (i.EQ.4) r2 = 5.0d1

            If (i.EQ.5) r1 = 5.0d1

            If (i.EQ.5) r2 = 1.0d2

            !write(*,*)i,' Lower,upper limits of radial grouping',r1,r2


            check1=.FALSE.

            If (rmax.LT.r2) check1=.TRUE.

            !write(*,*)i,' Does grid need reshaping?',r1,rmax,r2,check1


            If (check1) then


               do j=1,nc

                  do k=1,ng


                     r1=0.0d0

                     check2=.FALSE.

                     do l=1,nr+1

                        If (tfmask(i,k,l,j)) then

                           check2=.TRUE.

                           r1=tf(i,k,l,j,2)

                           exit

                        End If

                     end do

                     !write(*,*)i,j,k,' Inner boundary',r1,check2


                     r2=0.0d0

                     check4=.FALSE.

                     do l=nr+1,1,-1

                        If (tfmask(i,k,l,j)) then

                           check4=.TRUE.

                           r2=tf(i,k,l,j,2)

                           exit

                        End If

                     end do

                     !write(*,*)i,j,k,' Outer boundary',r2,check4

                     !write(50,*)i,j,k,r1,r2,check2,check4


                     check5=.FALSE.

                     check6=.FALSE.

                     check7=.FALSE.

                     If (check2.AND.check4) then

                        If ((r1.LT.rmax).AND.(r2.LT.rmax)) check5=.TRUE.

                        If ((r1.LT.rmax).AND.(r2.GT.rmax)) check6=.TRUE.

                        If ((r1.GT.rmax).AND.(r2.GT.rmax)) check7=.TRUE.

                     End If

                     !write(*,*)i,j,k,' Results of oundary tests',
!     &                                       check5,check6,check7

                     !write(51,*)i,j,k,r1,rmax,r2,check5,check6,check7


                     If (check6) then


                        do l=1,nr1

                           r(l)=tf(i,k,l,j,2)

                           m(l)=tf(i,k,l,j,3)

                           x(l)=tf(i,k,l,j,4)

                           y(l)=tf(i,k,l,j,5)

                           mask1a(l)=tfmask(i,k,l,j)

                        end do

                        !write(*,*)i,j,k,' Setup interpolation arrays'


                        scale=1.0d-2

                        call gridinit (z,r1,rmax,scale,nr1)

                        !write(*,*)i,j,k,'Initialised new radial grid'


                        do l=1,nr1

                           w(1)=0.0d0

                           w(2)=0.0d0

                           w(3)=0.0d0

                           mask1b(1)=.FALSE.

                           mask1b(2)=.FALSE.

                           mask1b(3)=.FALSE.

                           tf(i,k,l,j,2)=0.0d0

                           tf(i,k,l,j,3)=0.0d0

                           tf(i,k,l,j,4)=0.0d0

                           tf(i,k,l,j,5)=0.0d0

                           tfmask(i,k,l,j)=.FALSE.

                           call linterpolate (mask1b(1),w(1),z(l),
     &                                        mask1a,r,m,nr1)

                           call linterpolate (mask1b(2),w(2),z(l),
     &                                        mask1a,r,x,nr1)

                           call linterpolate (mask1b(3),w(3),z(l),
     &                                        mask1a,r,y,nr1)

                           If (mask1b(1).AND.mask1b(2)
     &                         .AND.mask1b(3)) then

                              tfmask(i,k,l,j)=.TRUE.

                              tf(i,k,l,j,2)=z(l)

                              tf(i,k,l,j,3)=w(1)

                              tf(i,k,l,j,4)=w(2)

                              tf(i,k,l,j,5)=w(3)

                           End If


                        end do


                     Else if (check7) then


                        do l=1,nr+1

                           tf(i,k,l,j,2)=0.0d0

                           tf(i,k,l,j,3)=0.0d0

                           tf(i,k,l,j,4)=0.0d0

                           tf(i,k,l,j,5)=0.0d0

                           tfmask(i,k,l,j)=.FALSE.

                        end do


                     End If


                  end do

               end do


            End If


            exit


         End If


      end do


      return

      end



      subroutine lcalc (checkout,g,nphot,dg,gmask,gpmask,tfmask,tf,
     &                         index,a,rtype,mtype,ngp,nc,ng,nr,np)

      implicit none

      logical checkout
      integer rtype,mtype
      double precision index,a
      integer ngp,nc,ng,nr,np
      double precision nphot(ngp,nc,ng)
      double precision g(ngp,nc,0:ng-1)
      double precision dg(ngp,nc,0:ng-1)
      double precision tf(ngp,ng+1,nr+1,nc,np)
      logical tfmask(ngp,ng+1,nr+1,nc)
      logical gmask(ngp,nc,0:ng-1)
      logical gpmask(ngp)

      double precision gval,xy(2,2,4),area,emm

      logical check1

      integer nvalid,counter
      integer i,j,k,l,m,n


      counter=0

      do i=1,ngp

         If (gpmask(i)) then

            do j=1,nc

               do k=1,ng


                  nphot(i,j,k) = 0.0d0

                  gmask(i,j,k-1)=.FALSE.

                  g(i,j,k-1) = tf(i,k,1,j,1)

                  If (k.EQ.ng) g(i,j,k-1)=tf(i,k+1,1,j,1)

                  dg(i,j,k-1) = tf(i,k+1,1,j,1)-tf(i,k,1,j,1)

                  gval = tf(i,k,1,j,1) + 0.5d0 * dg(i,j,k-1)

                  !write(*,*)i,j,k,' Redshift',gval


                  do l=1,nr


                     !write(*,*)i,j,k,l,' Starting iteration'


                     nvalid=0

                     do m=1,2

                        do n=1,2

                           xy(m,n,1)=tf(i,m+k-1,n+l-1,j,4)

                           xy(m,n,2)=tf(i,m+k-1,n+l-1,j,5)

                           xy(m,n,3)=tf(i,m+k-1,n+l-1,j,2)

                           xy(m,n,4)=tf(i,m+k-1,n+l-1,j,3)

                           If (tfmask(i,m+k-1,n+l-1,j)) nvalid=nvalid+1

                        end do

                     end do


                     !write(*,*)i,j,k,l,' Extracted box'


                     check1=.FALSE.

                     If (nvalid.EQ.4) check1=.TRUE.

                     !write(*,*)i,j,k,l,' Is box valid?',nvalid,check1


                     If (check1) then


                        call areacalc (area,xy)

                        !write(*,*)i,j,k,l,' Area of box',area


                        call emmcalc (emm,xy,index,a,rtype,mtype)

                        !write(*,*)i,j,k,l,' Emmisivity weighting',emm


                        counter=counter+1

                        gmask(i,j,k-1)=.TRUE.

                        nphot(i,j,k)=nphot(i,j,k)+
     &                            (gval*gval*gval*gval*area*emm)

                        !write(*,*)i,j,k,l,' Number of photons',nphot(i,j,k)

                        !write(16,*)i,j,k,l,gval,area,emm,nphot(i,j,k)


                     End If


                  end do

                  !write(17,*)i,j,k,g(i,j,k-1),nphot(i,j,k),dg(i,j,k-1)

                  !If (gmask(i,j,k-1)) !write(18,*)g(i,j,k-1),nphot(i,j,k)


               end do

               !write(18,*)'no no'

            end do

         End If

      end do


      checkout=.FALSE.

      If (counter.GT.0) checkout=.TRUE.

      !write(*,*)'Were valid points found in line profile?',checkout


      return

      end



      subroutine rebinp1 (nmask,n,e,de,emin,emax,gpmask,gmask,
     &                    nphot,dg,g,nb,ngp,ng,nc)

      implicit none

      integer nb,ngp,ng,nc
      double precision de(0:nb-1)
      double precision e(0:nb-1),n(nb)
      double precision nphot(ngp,nc,ng)
      double precision g(ngp,nc,0:ng-1)
      double precision dg(ngp,nc,0:ng-1)

      logical gpmask(ngp)
      logical nmask(0:nb-1)
      logical gmask(ngp,nc,0:ng-1)

      integer ng1
      parameter (ng1=100)

      double precision emin,emax,deval,norm
      double precision eval(ng1),nval(ng1),x,y

      logical check1
      logical mask1(ng1)

      integer i,j,k


c      emin=2.0d1

c      emax=-2.0d1

c      do i=1,ngp

c         If (gpmask(i)) then

c            do j=1,nc

c               emin=min(emin,g(i,j,0))

c               emax=max(emax,g(i,j,99))

c            end do

c         End If

c      end do

      !write(*,*)'Lower, upper limits of global grid',emin,emax


      e(0)=emin

      nmask(0)=.FALSE.

      deval=(emax-emin)/dble(nb-1)

      de(0)=deval

      do i=1,nb

         n(i)=0.0d0

         If (i.LT.nb) then

            de(i)=deval

            nmask(i)=.FALSE.

            e(i)=e(i-1)+deval

         End If

         !write(19,*)i,e(i-1),n(i),de(i-1),nmask(i-1)

      end do

      !write(*,*)'Initialised global grids'


      norm = 0.0d0

      do i = 1,ngp

         If (gpmask(i)) then

            do j=1,nc


               do k=1,ng

                  mask1(k)=gmask(i,j,k-1)

                  nval(k)=nphot(i,j,k)/dg(i,j,k-1)

                  eval(k)=g(i,j,k-1)+0.5d0*(dg(i,j,k-1))

                  !write(20,*)i,j,k,eval(k),nval(k),
!     &                       mask1(k),dg(i,j,k-1),gmask(i,j,k-1)

               end do

               !write(*,*)i,j,' Initialised interpolation grids'


               do k=1,nb

                 check1=.FALSE.

                 x=e(k-1)+0.5d0*deval

                 call linterpolate (check1,y,x,mask1,eval,nval,ng)

                 y = y * deval / x

                 If (check1) nmask(k)=.TRUE.

                 If (check1) n(k) = n(k) + y

                 If (check1) norm = norm + y

                 !write(21,*)i,j,k,x,y,n(k),check1

               end do


            end do

         End If

      end do

      !write(*,*)'Rebinned lines onto global grid'


      do i = 1,nb

         n(k) = n(k) / norm

         !write(22,*)i,e(i-1),n(i)

      end do

      !write(*,*)'Normalised line'


      return

      end



      subroutine lsmooth (nmask,n,e,de,gpmask,gmask,g,nb,ngp,ng,nc)

      implicit none

      integer nb,ngp,ng,nc
      double precision de(0:nb-1)
      double precision e(0:nb-1),n(nb)
      double precision g(ngp,nc,0:ng-1)

      logical gpmask(ngp)
      logical nmask(0:nb-1)
      logical gmask(ngp,nc,0:ng-1)

      double precision d1,d2,dnde,tol
      double precision gbound(ngp,nc,2)

      integer i,j,k


      tol=1.0d-1

      do i=1,ngp


         If (gpmask(i)) then


            do j=1,nc


               gbound(i,j,1)=g(i,j,0)

               gbound(i,j,2)=g(i,j,ng-1)

               !write(*,*)i,j,' Limits on line',gbound(i,j,1),
!     &                   gbound(i,j,2)


               do k=1,nb-4

                  If (e(k).LE.gbound(i,j,1).AND.
     &                e(k+1).GE.gbound(i,j,1)) then

                     d1=dabs(e(k)-gbound(i,j,1))

                     d2=dabs(e(k+1)-gbound(i,j,1))

                     !write(*,*)i,j,k,1,e(k-1),n(k),nmask(k-1)

                     !write(*,*)i,j,k,1,e(k),n(k+1),nmask(k),d1

                     !write(*,*)i,j,k,1,e(k+1),n(k+2),nmask(k+1),d2

                     If (d1.LT.d2) then

                        nmask(k)=.FALSE.

                        dnde=dabs((n(k+1)-n(k))/de(k-1))

                        If (dnde.LT.tol) nmask(k-1)=.FALSE.

                        !write(*,*)i,j,k,1,1,dnde,de(k-1)

                     Else

                        nmask(k)=.FALSE.

                        dnde=dabs((n(k+2)-n(k+1))/de(k+1))

                        If (dnde.LT.tol) nmask(k+1)=.FALSE.

                        !write(*,*)i,j,k,1,2,dnde,de(k+1),nmask(k+1)

                     End If

                  End If


                  If (e(k).LE.gbound(i,j,2).AND.
     &                e(k+1).GE.gbound(i,j,2)) then

                     d1=dabs(e(k)-gbound(i,j,2))

                     d2=dabs(e(k+1)-gbound(i,j,2))

                     !write(*,*)i,j,k,2,e(k-1),n(k),nmask(k-1)

                     !write(*,*)i,j,k,2,e(k),n(k+1),nmask(k),d1

                     !write(*,*)i,j,k,2,e(k+1),n(k+2),nmask(k+1),d2

                     If (d1.LT.d2) then

                        nmask(k)=.FALSE.

                        dnde=dabs((n(k+1)-n(k))/de(k-1))

                        If (dnde.LT.tol) nmask(k-1)=.FALSE.

                        !write(*,*)i,j,k,2,1,dnde,de(k-1)

                     Else

                        nmask(k)=.FALSE.

                        dnde=dabs((n(k+2)-n(k+1))/de(k+1))

                        If (dnde.LT.tol) nmask(k+1)=.FALSE.

                        !write(*,*)i,j,k,2,2,dnde,de(k+1),nmask(k+1)


                     End If

                  End If


               end do



            end do


         End If


      end do


      return

      end



      subroutine rebinp2 (ne,ec,ear,photar,nmask,n,e,de,nb)

      implicit none

      integer nb,ne
      real ear(0:ne),photar(ne)
      double precision ec,de(0:nb-1)
      double precision e(0:nb-1),n(nb)
      logical nmask(0:nb-1)

      double precision x,y,norm
      double precision eval(nb),nval(nb)

      logical check1
      logical mask1(nb)
      logical mask1a(ne)

      integer i


      do i=1,nb

         mask1(i)=nmask(i-1)

         nval(i)=n(i)/de(i-1)

         eval(i)=ec*(e(i-1)+0.5d0*(de(i-1)))

         !write(24,*)i,eval(i),nval(i),mask1(i)

      end do

      !write(*,*)'Initialised interpolation grids'


      norm=0.0d0

      do i=1,ne


         y = 0.0d0

         check1=.FALSE.

         photar(i)=0.0e0

         x=dble(ear(i-1)+0.5d0*(ear(i)-ear(i-1)))

         call linterpolate (check1,y,x,mask1,eval,nval,nb)

         y = y * (ear(i)-ear(i-1))

         If (check1) norm = norm + y

         If (check1) photar(i) = real(y)

         mask1a(i)=check1

         !write(25,*)i,x,y,photar(i),check1


      end do

      !write(*,*)'Rebinned line onto XSPEC grid',norm


      do i=1,ne

         If (mask1a(i)) photar(i) = photar(i) / real(norm)

c         write(26,*)i,ear(i-1),photar(i)

c         write(27,*)ear(i-1),photar(i)

      end do

      !write(*,*)'Normalised line'

c      write(26,*)

c      write(27,*)


      return

      end

      



      subroutine areacalc (area,xy)

      implicit none

      double precision area,xy(2,2,4)

      double precision dx,dy,ds
      double precision a,b,c,d
      double precision ang1,ang2,s
      double precision zero,one,tol

      logical mask1(5)
      logical check1,check2
      logical check4

      integer nvalid


      zero=0.0d0

      one=1.0d0

      tol=1.0d-10

      area=0.0d0

      !write(12,*)xy(1,1,1),xy(1,1,2)

      !write(12,*)xy(1,2,1),xy(1,2,2)

      !write(12,*)xy(2,2,1),xy(2,2,2)

      !write(12,*)xy(2,1,1),xy(2,1,2)

      !write(12,*)xy(1,1,1),xy(1,1,2)

      !write(12,*)'no no'


      nvalid=0

      mask1(1)=.FALSE.

      dx=dabs(xy(1,1,1)-xy(1,2,1))

      dy=dabs(xy(1,1,2)-xy(1,2,2))

      a=(dx*dx)+(dy*dy)

      If (a.GT.tol) then

         nvalid=nvalid+1

         mask1(1)=.TRUE.

         a=dsqrt(a)

      End If

      !write(*,*)'Length of side a',a,dx,dy


      mask1(2)=.FALSE.

      dx=dabs(xy(1,1,1)-xy(2,1,1))

      dy=dabs(xy(1,1,2)-xy(2,1,2))

      b=(dx*dx)+(dy*dy)

      If (b.GT.tol) then

         nvalid=nvalid+1

         mask1(2)=.TRUE.

         b=dsqrt(b)

      End If

      !write(*,*)'Length of side b',b,dx,dy


      mask1(3)=.FALSE.

      dx=dabs(xy(2,1,1)-xy(2,2,1))

      dy=dabs(xy(2,1,2)-xy(2,2,2))

      c=(dx*dx)+(dy*dy)

      If (c.GT.tol) then

         nvalid=nvalid+1

         mask1(3)=.TRUE.

         c=dsqrt(c)

      End If

      !write(*,*)'Length of side c',c,dx,dy


      mask1(4)=.FALSE.

      dx=dabs(xy(1,2,1)-xy(2,2,1))

      dy=dabs(xy(1,2,2)-xy(2,2,2))

      d=(dx*dx)+(dy*dy)

      If (d.GT.tol) then

         nvalid=nvalid+1

         mask1(4)=.TRUE.

         d=dsqrt(d)

      End If

      !write(*,*)'Length of side d',d,dx,dy


      mask1(5)=.FALSE.

      dx=dabs(xy(1,1,1)-xy(2,2,1))

      dy=dabs(xy(1,1,2)-xy(2,2,2))

      ds=(dx*dx)+(dy*dy)

      If (ds.GT.tol) then

         nvalid=nvalid+1

         mask1(5)=.TRUE.

         ds=dsqrt(ds)

      End If

      !write(*,*)'Length of base of triangles',ds


      check1=.FALSE.

      If (nvalid.EQ.5) check1=.TRUE.

      !write(*,*)'Are all lengths valid?',nvalid,check1

      !write(14,*)a,b,c,d,ds,nvalid,check1


      If (check1) then


         check2=.FALSE.

         ang1=(a*a)+(d*d)-(ds*ds)

         ang1=ang1/((a*d)+(a*d))

         If (dabs(ang1).LE.one) check2=.TRUE.

         If (check2) ang1=dacos(ang1)

         !write(*,*)i,j,'Value of first angle',ang1


         check4=.FALSE.

         ang2=(b*b)+(c*c)-(ds*ds)

         ang2=ang2/((b*c)+(b*c))

         If (dabs(ang2).LE.one) check4=.TRUE.

         If (check4) ang2=dacos(ang2)

         !write(*,*)i,j,'Value of second angle',ang2



         If (check2.AND.check4) then


            s=a+b+c+d

            s=0.5d0*s

            !write(*,*)i,j,'Length of semiperimeter',s

            !write(19,*)ang1,ang2,s


            area=dcos(0.5d0*(ang1+ang2))

            area=area*area

            area=a*b*c*d*area

            area=((s-a)*(s-b)*(s-c)*(s-d))-area

            area=dsqrt(area)

            !write(*,*)i,j,'Area of trapezoid',area

            !write(20,*)area


          End If


      Else


        nvalid=0

c        do i=1,4

c          If (mask1(i)) nvalid=nvalid+1

c        end do

        !write(*,*)'Number of valid edges',nvalid


        check2=.FALSE.

        If (nvalid.EQ.3) check2=.TRUE.

        !write(*,*)'Is shape triangular?',check2


        If (check2) then


          If (.NOT.(mask1(1))) then

             a=b

             b=c

             c=d

          Else if (.NOT.(mask1(2))) then

             b=c

             c=d

          Else if (.NOT.(mask1(3))) then

             c=d

          End If

          !write(*,*)'Sides of triangle',a,b,c


          area=(a+b+c)*(b+c-a)

          area=area*(c+a-b)*(a+b-c)

          area=dsqrt(area)

          area=0.25d0*area

          !write(*,*)'Area of triangles',area


       End If




      End If


      return

      end



      subroutine emmcalc (emm,xy,id,a,rt,mt)

      implicit none

      integer rt,mt
      double precision emm,xy(2,2,4),a,id

      double precision re,me,fztibc

      integer i,j


      re=0.0d0

      me=0.0d0

      do i=1,2

         do j=1,2

            re=re+xy(i,j,3)

            me=me+xy(i,j,4)

         end do

      end do

      re = 0.25d0 * re

      me = 0.25d0 * me

      !write(*,*)'Averaged values of re, me',re,me


      emm=1.0d0

      If (rt.EQ.1) then

        emm = emm * (re**(-id))

      Else if (rt.EQ.2) then

        emm = emm * fztibc(re,a)

      End If


      If (mt.EQ.1) then

        emm = emm / me

      Else if (mt.EQ.2) then

        emm = emm * ((1.0d0+(2.06d0*me))/2.06d0)

      Else if (mt.EQ.3) then

        emm = emm * (me*dlog(1.0d0+(1.0d0/me)))

      End If

      !write(*,*)'Emissivity',emm

      !write(15,*)re,me,id,emm


      return

      end



      subroutine pathinit (path,rmin,rmax,a,obs)

      double precision rmin,rmax,a,obs

      character*256 path,fgmodf

      character*21 stringp1
      character*36 stringp2
      character*4 s1
      character*3 s2
      character*9 fmt

      integer x(4),zero
      integer pathlength
      integer lenact

      external fgmodf,lenact


      zero=0

      stringp1='0000000000000.FITS.gz'

      !write(*,*)'Input variables',a,obs,rmin,rmax

      !write(*,*)'Initialised stringp1 ',stringp1


      x(1)=nint(a*1.0d3)

      x(2)=nint(obs)

      x(3)=nint(rmin)

      x(4)=nint(rmax)

      !write(*,*)'Setup values',x


      fmt='(I4)'

      s1 = ' '

      write(s1,fmt) x(1)

      stringp1(1:4) = s1

      !write(*,*)'First component of file name ',stringp1


      fmt='(I3)'

      s2 = ' '

      write(s2,fmt) x(2)

      stringp1(5:7) = s2

      !write(*,*)'Second component of file name ',stringp1


      fmt='(I3)'

      s2 = ' '

      write(s2,fmt) x(3)

      stringp1(8:10) = s2

      !write(*,*)'Third component of file name ',stringp1


      fmt='(I3)'

      s2 = ' '

      write(s2,fmt) x(4)

      stringp1(11:13) = s2

      !write(*,*)'Fourth component of file name ',stringp1


      fmt='(I1)'

      If (x(1).LT.1000) write(stringp1(1:1),fmt) zero

      If (x(1).LT.100) write(stringp1(2:2),fmt) zero

      If (x(1).LT.10) write(stringp1(3:3),fmt) zero

      If (x(2).LT.100) write(stringp1(5:5),fmt) zero

      If (x(2).LT.10) write(stringp1(6:6),fmt) zero

      If (x(3).LT.100) write(stringp1(8:8),fmt) zero

      If (x(3).LT.10) write(stringp1(9:9),fmt) zero

      If (x(4).LT.100) write(stringp1(11:11),fmt) zero

      If (x(4).LT.10) write(stringp1(12:12),fmt) zero

      !write(*,*)'Correct filename ',stringp1


      stringp2(1:14)='XSKDLINE/0000/'

      stringp2(10:13) = s1

      If (x(1).LT.1000) write(stringp2(10:10),fmt) zero

      If (x(1).LT.100) write(stringp2(11:11),fmt) zero

      If (x(1).LT.10) write(stringp2(12:12),fmt) zero

      stringp2(15:36)=stringp1

c      write(*,*)'Path to file ',stringp2,stringp1


      path = fgmodf()

c      write(*,*)'Path to Data directory',datadir


      pathlength = lenact(path)

c      write(*,*)'Length of path'


      path(pathlength+1:pathlength+37) = stringp2

c      write(*,*)'Full path to file',path


      return

      end



      subroutine tfopen (checkout,mask,g,r,m,x,y,path,ng,nr,nc)

      implicit none

      integer ng,nr,nc
      character*256 path
      double precision g(ng+1,nr+1,nc)
      double precision r(ng+1,nr+1,nc)
      double precision m(ng+1,nr+1,nc)
      double precision x(ng+1,nr+1,nc)
      double precision y(ng+1,nr+1,nc)
      logical checkout,mask(ng+1,nr+1,nc)

      integer n
      parameter (n=20402)

      double precision nullval
      double precision a(n),b(n),c(n),d(n),e(n)

      logical m1(n),check1

      character errtext*30,errmessage*80

      integer unit,status,readwrite
      integer blocksize,hdutype
      integer frow,felem,colnum
      integer i,j,k,l


      status=0

      call ftgiou(unit,status)

c      write(*,*)'Got unit number',unit,status


      readwrite=0

      call ftopen(unit,path,readwrite,blocksize,status)

c      write(*,*)'Opened file, status',status


      call ftmahd(unit,2,hdutype,status)

c      write(*,*)'Moved to binary table',hdutype,status


      frow=1

      felem=1

      colnum=1

      nullval=0.0d0

      check1=.FALSE.

      call ftgcvd (unit,colnum,frow,felem,n,nullval,a,check1,status)

c      write(*,*)'Read in redshifts, result',check1,status


      colnum=2

      check1=.FALSE.

      call ftgcvd (unit,colnum,frow,felem,n,nullval,b,check1,status)

c      write(*,*)'Read in radii, result',check1,status


      colnum=3

      check1=.FALSE.

      call ftgcvd (unit,colnum,frow,felem,n,nullval,c,check1,status)

c      write(*,*)'Read in angles, result',check1,status


      colnum=4

      check1=.FALSE.

      call ftgcvd (unit,colnum,frow,felem,n,nullval,d,check1,status)

c      write(*,*)'Read in alpha coords, result',check1,status


      colnum=5

      check1=.FALSE.

      call ftgcvd (unit,colnum,frow,felem,n,nullval,e,check1,status)

c      write(*,*)'Read in beta coords, result',check1,status


      colnum=6

      check1=.FALSE.

      call ftgcl (unit,colnum,frow,felem,n,m1,check1,status)

c      write(*,*)'Read in mask, result',check1,status


      call ftclos(unit, status)

      call ftfiou(unit, status)

      If (status.GT.0) then


         call ftgerr(status,errtext)

         write(*,*)'FITSIO Error status =',status,': ',errtext


         call ftgmsg(errmessage)

         do while (errmessage.NE. ' ')

            write(*,*)errmessage

            call ftgmsg(errmessage)

         end do


      End If

c      write(*,*)'Status of FITS read',status


      checkout=.FALSE.

      If (status.EQ.0) checkout=.TRUE.

c      write(*,*)'Was FITS read succesfull?',checkout


      If (checkout) then


         i=1

         do j=1,ng+1

            do k=1,nr+1

               do l=1,nc

                  g(j,k,l)=a(i)

                  r(j,k,l)=b(i)

                  m(j,k,l)=c(i)

                  x(j,k,l)=d(i)

                  y(j,k,l)=e(i)

                  mask(j,k,l)=m1(i)

                  i=i+1

               end do

            end do

         end do

c         write(*,*)'Copied TF to output arrays'


      End If


      return

      end



      subroutine gridinit (x,xmin,xmax,scale,n)

      implicit none

      integer n
      double precision x(n),xmin,xmax,scale

      double precision step

      integer i


      x(1)=xmin

      x(n)=xmax

      step=(xmax-xmin)*scale

      do i=2,n-1

         x(i)=xmin+(dble(i-1)*step)

      end do


      return

      end




      double precision function rms (a)

      implicit none

      double precision a

      double precision a2,z1,z2
      double precision one,three,third


      one=1.0d0

      three=one+one+one

      third=one/three

      a2=a*a

      !write(*,*)'Setup parameters',one,three,third,a,a2


      z1=((one+a)**third)+((one-a)**third)

      z1=((one-a2)**third)*z1

      z1=z1+one

      !write(*,*)'Value of Z1',z1


      z2=a2+a2+a2+(z1*z1)

      z2=dsqrt(z2)

      !write(*,*)'Value of Z2',z2


      rms=(three-z1)*(three+z1+z2+z2)

      rms=dsqrt(rms)

      rms=three+z2-rms

      !write(*,*)'Location of rms',rms


      return

      end



      subroutine linterpolate (checkout,ytarg,xtarg,xymask,xin,yin,n)

      implicit none

      integer n
      double precision xtarg,ytarg
      double precision xin(n),yin(n)
      logical xymask(n),checkout

      double precision x(n),y(n),frac

      integer nvalid
      integer i,j


      i=1

      nvalid=0

      do j=1,n

        x(j)=0.0d0

        y(j)=0.0d0

        If (xymask(j)) then

          x(i)=xin(j)

          y(i)=yin(j)

          nvalid=nvalid+1

          i=i+1

        End If

      end do

      !write(*,*)'Sorted points into ascending order in x'

      !write(*,*)'Number of valid points',nvalid


      ytarg=0.0d0

      do i=1,nvalid-1


        If ((x(i).LE.xtarg).AND.(x(i+1).GE.xtarg)) then

          checkout=.TRUE.

          frac=(xtarg-x(i))/(x(i+1)-x(i))

          ytarg=y(i)+frac*(y(i+1)-y(i))

          exit

        End If


      end do

      !write(*,*)'Was valid point found?',xtarg,ytarg,checkout


      return

      end


      double precision function fztibc (r,a)

      implicit none

      double precision r,a

      double precision rms,x,pi
      double precision x0,x1,x2,x3
      double precision dm1,dm2,dm3
      double precision dm4,dm5,dm6
      double precision third


      pi = dacos(-1.0d0)

      third = 1.0d0/3.0d0


      x = dsqrt(r)

      x0 = dsqrt(rms(a))

      x1 = 2.0d0*dcos(third*(dacos(a) - pi))

      x2 = 2.0d0*dcos(third*(dacos(a) + pi))

      x3 = -2.0d0*dcos(third*dacos(a))


      dm1 = 0.5d0*3.0d0*a*dlog(x/x0)


      dm2 = dlog((x-x1)/(x0-x1))

      dm2 = 3.0d0*(x1-a)*(x1-a)*dm2

      dm2 = dm2 / (x1*(x1-x2)*(x1-x3))


      dm3 = dlog((x-x2)/(x0-x2))

      dm3 = 3.0d0*(x2-a)*(x2-a)*dm3

      dm3 = dm3 / (x2*(x2-x1)*(x2-x3))


      dm4 = dlog((x-x3)/(x0-x3))

      dm4 = 3.0d0*(x3-a)*(x3-a)*dm4

      dm4 = dm4 / (x3*(x3-x1)*(x3-x2))


      dm5 = x - x0 - dm1 - dm2 - dm3 - dm4


      dm6 = (x*x*x) - x - x - x + a + a

      dm6 = x * x * dm6

      dm6 = 0.5d0 * 3.0d0 *dm5 / dm6


      fztibc = dm6 / r


      end function
