c-----------------------------------------------------------------------
c metric.f
c  Written by Mijan Huq
c   Center for Relativity,
c   University of Texas at Austin.
c
c	This file contains routines that calculate the components of the
c	Novikov metric
c-----------------------------------------------------------------------
c Routine to calculate metric coefficients for Novikov data.
c Calls the routine NewtSlv to solve for eta on a tau-Rstar
c grid. Then this routine goes about generating r and t. Difference
c r and t to generate metric components.
c	Modification MFH
c		o Replaced finite differencing of r and t by analytic expressions
c		  for gRR, gtt, KRR and Ktt. t:=theta here. Cleaned up code...
c
c
		subroutine NovMetric(tau, Rstar, eta, r, t, gRR, gtt,kRR,ktt,
     .	detadtau,detadRstar,
     .	drdtau,drdRstar,dtdtau,dtdRstar, M, stp, itmax,NT,NR)
			implicit none

			integer	NT, NR, itmax

			real	M, stp
			real	tau(NT), Rstar(NR), eta(NT,NR),r(NT,NR), t(NT,NR), 
     .		drdtau(NT,NR), drdRstar(NT,NR), detadRstar(NT,NR),
     .		dtdtau(NT,NR),dtdRstar(NT,NR), detadtau(NT,NR),
     .		gRR(NT,NR), gtt(NT,NR), kRR(NT,NR), ktt(NT,NR)

c		Local Variables

			integer	i, j, k

			real		Fact
			real		dtau, dRstar, d2R
		
			dtau = tau(2) - tau(1)

c		First thing to be done is to solve for eta

			write(0,*)'Solving for eta'
			call NewtSlv(Rstar,eta,tau,M,stp,itmax,NT,NR)

c		Here is another option for generating eta from a table.
c			call  GetEta(Rstar,eta,tau,M,NT,NR)

c		Now generate r and t
		
			write(0,*)'generating r and t'	
			call GetRT(tau,Rstar,eta,r,t,M,NT,NR)

c		Consistency check from MTW 31.12b

			write(0,*)'Consistency check'
			call checkme(tau,Rstar,r,M,NT,NR)

c		Calculate partials

			write(0,*)'Getting partials'
			call GetPartials(tau,Rstar,r,M,drdtau,drdRstar,NT,NR)

c		Calculate metric components gRR and gtt for each grid point.

			do i = 1, NT
				do j = 1, NR

					gRR(i,j) = (2.0e0 * r(i,j) + 3.0e0 * tau(i)*
     .			sqrt(2.0e0*M/r(i,j) - 1.0e0/(Rstar(j)**2 + 1.0e0)))**2/
     .			(Rstar(j)**2 + 1.0e0)

					gtt(i,j) = r(i,j)**2

				end do
			end do

c		Form ktt 


c		
			write(0,*)'Calculating kRR and ktt'
			do i = 1, NT
				do j = 1, NR

					Fact = sqrt(2.0e0 * M / r(i,j) - 
     .						1.0e0/(Rstar(j)**2 + 1.0e0))
					
c					if(tau(j) .gt. 1.0e-13)then
c					kRR(i,j) = 1.0e0 / (Rstar(j)**2 + 1.0e0) *
c     .				(2.0e0*r(i,j) - 3.0e0 *tau(i)*Fact)*
c     .				( Fact + 3.0e0*tau(i)*M/
c     .				r(i,j)**2 / Fact)
c					else
c					kRR(i,j) = 0.0e0
c					end if

					kRR(i,j) = -1.0e0/(Rstar(j)**2 + 1.0e0) *(2.0e0 *r(i,j)
     .					+ 3.0e0*tau(i)*Fact)*(Fact + 3.0e0*M*tau(i)/
     .					r(i,j)**2)

					ktt(i,j) = r(i,j)*Fact
					write(99,*)tau(i), Rstar(j), ktt(i,j)
					write(98,*)tau(i), Rstar(j), kRR(i,j)

				end do
				write(98,*)
				write(99,*)
			end do	

1011		continue
c			call Consist(tau,Rstar,eta,detadtau,r,M,NT,NR)

			return
		end
c-----------------------------------------------------------------------
c	Routine to generate r and t given eta on a tau-Rstar grid.
c
		subroutine GetRT(tau,Rstar,eta,r,t,M,NT,NR)
			implicit none

			integer	NT, NR

			real	tau(NT), Rstar(NR),
     .			eta(NT,NR), r(NT,NR), t(NT,NR),
     .			M

c		Local Variables

			integer	i, j, k

c		Loop over points in tau-Rstar grid


			do i = 1, NT
				do j = 1, NR

					r(i,j) = (Rstar(j)**2 + 1.0e0)*(1.0e0 + cos(eta(i,j)))*M

c					t(i,j) = 2.0e0*M*(log(abs(Rstar(j)+tan(eta(i,j)*0.5e0))/
c     .		abs(Rstar(j) - tan(eta(i,j)*0.5e0)))) +
c     .		2.0e0*M*Rstar(j)*(eta(i,j) + 0.5e0*(Rstar(j)**2 + 1.0e0)*
c     .		(eta(i,j) + sin(eta(i,j))))

				end do
			end do
			return
		end	
c-----------------------------------------------------------------------
c	Routine to difference r and t wrt tau and Rstar.
			subroutine DiffRT(tau,Rstar,r,t,
     .					drdtau,drdRstar,dtdtau,dtdRstar,NT,NR)
				implicit none

				integer	NT, NR

				real	tau(NT), Rstar(NR), r(NT,NR), t(NT,NR),
     .      		drdtau(NT,NR), drdRstar(NT,NR),
     .				dtdtau(NT,NR),dtdRstar(NT,NR)

c			Local Variables

				integer	i, j, k
	
				real	dtau, dRstar, divtau, divRstar

				dtau = tau(2) - tau(1)

				dRstar = Rstar(2) - Rstar(1)


c			Difference  r first

				call DiffRR(r,drdRstar,dRstar,NT,NR)

				call DiffTT(r,drdtau, dtau, NT,NR)

c			Difference t

				call DiffRR(t,dtdRstar,dRstar,NT,NR)

				call DiffTT(t,dtdtau, dtau, NT,NR)
				
			return
		end
c-----------------------------------------------------------------------
c	
c	The following routine carries out a consistency check on the
c	grid functions.  Using maple the terms were cast into fortran code.
c	See notes page 10.
c
		subroutine Consist(tau,Rstar,eta,detadtau,r,M,NT,NR)
			implicit none

			integer	NT, NR

			real	eta(NT,NR), detadtau(NT,NR), r(NT,NR),
     .			Rstar(NR), tau(NT),
     .			M

c		Local Variables

			integer	i,j,k

			real	term1a, term2a, term1b, term2b

			do i = 1, NT
				do j = 1, NR

c					drdeta*drdR / A term

      term1a = -2*M**2*sin(eta(i,j))*Rstar(j)*(Rstar(j)**4+2*Rstar(j)**4
     #		*cos(eta(i,j))+Rstar(j)**4*cos(eta(i,j))**2+2*Rstar(j)
     #**2+4*Rstar(j)**2*cos(eta(i,j))+2*Rstar(j)**2*cos(eta(i,j))**2
     #		+1+2*cos(eta(i,j))+cos(eta(i,j))**2)/(
     #	Rstar(j)**2+Rstar(j)**2*cos(eta(i,j))-1+cos(eta(i,j)))

c					dtdeta*dtdR*A term.

      term2a = (1-2/(Rstar(j)**2+1)/(1+cos(eta(i,j))))*(2*M*
     #		sign(1.E0,(Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-t
     #an(eta(i,j)/2)))*((1.E0/2.E0+tan(eta(i,j)/2)**2/2)/(Rstar(j)
     #		-tan(eta(i,j)/2))-(Rstar(j)+tan(eta(i,j)
     #/2))/(Rstar(j)-tan(eta(i,j)/2))**2*(-1.E0/2.E0-tan(eta(i,j)/2
     #		)**2/2))/abs((Rstar(j)+tan(et
     #	a(i,j)/2))/(Rstar(j)-tan(eta(i,j)/2)))+2*M*Rstar(j)*(1+
     #		(Rstar(j)**2+1)*(1+cos(eta(i,j)))/2))*(2*M*sign(
     #1.E0,(Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-tan(eta(i,j)/2)))*
     #		(1/(Rstar(j)-tan(eta(i,j)/2))-(Rstar(j)+tan(eta(i,j)/2
     #))/(Rstar(j)-tan(eta(i,j)/2))**2)/abs((Rstar(j)+tan(eta(i,j)
     #		/2))/(Rstar(j)-tan(eta(i,j)/2)))+2*M*(eta(i,j)+
     #(Rstar(j)**2+1)*(eta(i,j)+sin(eta(i,j)))/2)+2*M*Rstar(j)**2*
     #		(eta(i,j)+sin(eta(i,j))))

c					form drdeta*drdR / A  - dtdeta*dtdR*A 
c					Which should be zero to get to the Novikov metric

					write(1,*)tau(i),Rstar(j), term1a - term2a

c					A*dtdeta**2 term

      term1b = (2*M*sign(1.E0,(Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-
     #	tan(eta(i,j)/2)))*((1.E0/2.E0+tan(
     #	eta(i,j)/2)**2/2)/(Rstar(j)-tan(eta(i,j)/2))-(Rstar(j)+
     #	tan(eta(i,j)/2))/(Rstar(j)-tan(eta(i,j)/2))**2*(-1.E
     #0/2.E0-tan(eta(i,j)/2)**2/2))/abs((Rstar(j)+tan(eta(i,j)/2))/
     #	(Rstar(j)-tan(eta(i,j)/2)))+2*M*Rstar(j)*
     #(1+(Rstar(j)**2+1)*(1+cos(eta(i,j)))/2))**2*(1-2/(Rstar(j)**2
     #	+1)/(1+cos(eta(i,j))))

c					drdeta**2 /A
      term2b= -M**2*(-Rstar(j)**6+Rstar(j)**6*cos(eta(i,j))**2-
     #		cos(eta(i,j))*Rstar(j)**6+cos(eta(i,j))**3*Rstar(j)**6-
     #	3*Rstar(j)**4+3*Rstar(j)**4*cos(eta(i,j))**2-3*Rstar(j)**4
     #	*cos(eta(i,j))+3*cos(eta(i,j))**3*Rstar(j)**4-3*Rstar(j)**
     #2+3*Rstar(j)**2*cos(eta(i,j))**2-3*Rstar(j)**2*cos(eta(i,j))+
     #		3*cos(eta(i,j))**3*Rstar(j)**2-1+cos(eta(i,j))
     #**2-cos(eta(i,j))+cos(eta(i,j))**3)/(Rstar(j)**2+Rstar(j)**2
     #		*cos(eta(i,j))-1+cos(eta(i,j)))

c				Form detadtau^2 *(term1b - term 2b)

					write(2,*)tau(i),Rstar(j),detadtau(i,j)*(term1b -
     .					term2b) - 1.0e0

					write(3,*)tau(i), Rstar(j),
     .			4*M**2*Rstar(j)**2*(1+cos(eta(i,j)))**2/(1-2/
     .			(Rstar(j)**2+1)/(1+cos(eta(i,j))))
     .	- 		(2*M*sign(1.E0,(Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-
     #	tan(eta(i,j)/2)))*(1/(Rstar(j)-tan(eta(i,j)/2)
     #)-(Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-tan(eta(i,j)/2))**2)/
     #	abs((Rstar(j)+tan(eta(i,j)/2))/(Rstar(j)-tan(eta(i,j)/
     #2)))+2*M*(eta(i,j)+(Rstar(j)**2+1)*(eta(i,j)+sin(eta(i,j)))/2)
     #	+2*M*Rstar(j)**2*(eta(i,j)+sin(eta(i,j))))*
     #*2*(1-2/(Rstar(j)**2+1)/(1+cos(eta(i,j))))


				end do
				write(1,*)
				write(2,*)
			end do

			return
		end
c-----------------------------------------------------------------------
		subroutine GetPartials(tau,Rstar,r,M,drdtau,drdRstar,NT,NR)
			implicit none

			integer	NT, NR

			real	tau(NT), Rstar(NR), r(NT,NR), 
     .			drdtau(NT,NR), drdRstar(NT,NR),
     .			M

c		Local Variables

			integer	i, j, k

			real	A

			do i = 1, NT
				do j = 1, NR

					drdtau(i,j) = -sqrt(-1.0e0/(Rstar(j)**2 + 1.0e0) + 
     .								2.0e0*M/r(i,j) )/ M

					drdRstar(i,j) = Rstar(j)/(Rstar(j)**2 + 1.0e0) *
     .			(2.0e0 * r(i,j)  + 3.0e0 * tau(i)*sqrt(
     .			-1.0e0/(Rstar(j)**2 + 1.0e0) + 2.0e0*M/r(i,j) ))

c					write(20,*)r(i,j),drdtau(i,j),drdRstar(i,j)
c					write(21,*)-1.0e0/(Rstar(j)**2 + 1.0e0) + 2.0e0*M/r(i,j) 

				end do
c				write(20,*)
			end do

			return
		end
c-----------------------------------------------------------------------
		subroutine checkme(tau,Rstar,r,M,NT,NR)
			implicit none

			integer	NT, NR

			real	tau(NT), Rstar(NR), r(NT,NR), M

c		Local variables

			integer	i,j,k

			real		F	, Fact1, Fact2, Fact3

			do i = 1, NT
				do j = 1, NR

					Fact1 = Rstar(j)**2 + 1.0e0
					Fact2 = r(i,j)/M
					Fact3 = Fact2 / Fact1 - 1.0e0

c		NOTE patched negative square root with abs statment.
c		Fact2 - (Fact2**2) / Fact1 was on the order of -1.0e-7 on the
c		level of roundoff error. Got negative roundoff error and so sqrt
c		gave nan's. -MFH

					if((Fact2/Fact1 - 1.0e0) .gt. 1.0e0)then
c						write(0,*)sqrt(Fact2/Fact1),' at ',tau(i),Rstar(j)
						Fact3 = 1.0e0
c						stop
					end if

					F = 0.5e0 * Fact1**1.5 * acos(Fact3) +
     .				Fact1 * sqrt(abs(0.25e0 * Fact2**2 / 
     .				Fact1 - 0.5e0*Fact2)) - tau(i)/2.0e0/M

					write(3,*)tau(i),Rstar(j),F


				end do
				write(3,*)
			end do

			return
		end
c-----------------------------------------------------------------------