c-----------------------------------------------------------------------
c newton.f 
c Created Fri Apr 14 15:55:13 CDT 1995 MFH
c
c	In this file there are two sets of routines for generating eta as
c	a function of tau and Rstar. The first implements a newton's method
c	and the second evaluates sin(eta) + eta = f and stores values of
c	eta and corresponding values of f in another array. Then an interpolation
c	can be carried out in f to get eta. This should work since f is a 
c	single valued function of eta.
c-----------------------------------------------------------------------
c	Modification Thu Apr 20 14:05:50 CDT 1995
c	Complete modification of code to solve different problem.
c	Instead of solving previous form of equation to get r as a function
c	of tau and Rstar we solve for eta as a function of tau and Rstar.
c	Then we go about constructing r and t. 
c
c
c	stp:= Stopping Criterion

		subroutine NewtSlv(Rstar,eta,tau,M,stp,itmax,NT,NR)
			implicit none

			integer NT, NR, itmax

			real	Rstar(NR), eta(NT,NR), tau(NT)
			real	M, stp

c	Local Variables

			integer i,j,k, ll

			real	F, dF, deltaeta, eta0

c	Functions called
			real ff, Dff, num

c	Set up initial eta0

			eta0 = 0.0e0



c		For each tau time step
			do i = 1, NT
c			And for each point in Rstar
				do j = 1, NR

					do k = 1, itmax
c					Use r0 as a guess	
						F = ff(eta0, Rstar(j), tau(i),M)
						dF = Dff(eta0, Rstar(j), tau(i),M)

						deltaeta = -F / dF

						eta0 = eta0 + deltaeta

c						write(0,*)eta0,F,dF,deltaeta



						if(abs(F) .lt. stp)then
							goto 201
						end if
					end do
					if(abs(F) .gt. stp)then
						write(0,*)'Newton method did not converge to',stp
						write(0,*)'NT = ',i, 'NR = ',j
						write(0,*)'F = ',F
						write(0,*)'STOPPING'
						stop
					end if
201		continue
c						write(0,*)tau(i),Rstar(j),F,dF,eta0

					eta(i,j) = eta0

c					Use eta0 as initial guess on next iteration
				end do
c					At this point set eta back to first grid point in R
c					since that is where we will start with the j loop
				eta0 = eta(i,1)
			end do

c			write(0,*)'Dumping eta'

c			do i = 1, NT
c				do j = 1, NR
c					write(80,*)tau(i),Rstar(j),eta(i,j)
c     				write(81,*)tau(i),Rstar(j),ff(eta(i,j),Rstar(j),tau(i),
c     .							M),
c     .				Dff(eta(i,j),Rstar(j),tau(i),M)
c				end do
c					write(80,*)
c					write(81,*)
c			end do


			return
		end
c-----------------------------------------------------------------------
c		Routines ff and Dff checked with maple-MFH
c			everything kosher so far

		real function ff(eta,Rstar,tau,M)
			implicit none

			real eta, Rstar, tau, M
			
				ff = eta + sin(eta) - tau / M / (Rstar**2 + 1)**1.5
			return
		end 
		
c-----------------------------------------------------------------------

		real function Dff(eta,Rstar,tau,M)
			implicit none

			real eta, Rstar, tau, M, rr

				Dff = 1.0e0 + cos(eta)

			return
		end 
		
c-----------------------------------------------------------------------
c		Initialization of coordinates 
c
		subroutine InitCrd(Rstar, tau,Rmin, Rmax, tmin, tmax,NR,NT)
			implicit none

			integer NR, NT

			real	Rstar(NR), tau(NT)

			real	Rmin,Rmax,tmin,tmax

c		Local Variables

			integer	i,j,k

			real	dRstar, dtau

c		Set up grid spacings
			dRstar = (Rmax - Rmin) / float(NR-1)
			dtau   = (tmax - tmin) / float(NT-1)

			write(0,*)'Grid Spacings',NR,NT
			write(0,*)'Space : ',dRstar
			write(0,*)'Time : ',dtau

			write(0,*)'Rmax = ',Rmax
			write(0,*)'Rmin = ',Rmin
			write(0,*)'tmax = ',tmax
			write(0,*)'tmin = ',tmin

c		Set up Rstar coordinates
			do i = 1 , NR
				Rstar(i) = Rmin + float(i-1)*dRstar
			end do

c		Set up tau coordinate
			do i = 1, NT
				tau(i) = tmin + float(i-1) * dtau
			end do

			return
		end
c=======================================================================
c	Begining of second set of routines to get eta.
c
c		Routine to generate table of values for sin x + x = f
c
		subroutine mktable(eta,rhs,N)
			implicit none

			integer	N
	
			real	eta(N), rhs(N)

c		Local Variables

			integer	i, j, k

			real	deta

			deta = 4.0e0 * atan(1.0e0) / float(N-1)

			do i = 1, N
				eta(i) = float(i-1)*deta 
				rhs(i) = eta(i)  + sin(eta(i))
			end do

			return
		end
c-----------------------------------------------------------------------
c
		subroutine intptable(eta,rhs,findme,thisrhs,N)
			implicit none
			
			integer	N

			real	eta(N), rhs(N), findme, thisrhs

c		Local variable 

			integer	it, i, j

			real	xx(2), ff(2)

c		Find it in the non-uniform grid in rhs

			do i = 1, N-1
c				if(rhs(i+1) .ge. thisrhs .and. rhs(i) .le. thisrhs)then
				if(rhs(i+1) .ge. thisrhs)then
					it = i
					goto 1287
				end if
			end do

1287		continue


			ff(1) = eta(it)
			ff(2) = eta(it + 1)

			xx(1) = rhs(it)
			xx(2) = rhs(it + 1)

c			write(0,*)thisrhs - xx(1)
c			write(0,*)thisrhs - xx(2)

c			write(0,*)xx(1), thisrhs, xx(2)
c			write(0,*)ff(1), ff(2)

			findme = (thisrhs - xx(2))/(xx(1) - xx(2)) * ff(1) +
     .				(thisrhs - xx(1))/(xx(2) - xx(1)) * ff(2)

			return
		end
c			
c-----------------------------------------------------------------------
c
		subroutine GetEta(Rstar,eta,tau,M,NT,NR)
			implicit none

			integer	NT, NR

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

c		Local Variables

			integer	i, j, k, NETA

			real		leta(100), lrhs(100), teta, trhs

			NETA = 100

c		Generate eta table
			call mktable(leta,lrhs,NETA)

			do i = 1, NT
				do j = 1, NR
					trhs = tau(i) / M / (Rstar(j)**2 + 1.0e0)**1.5	
c		Find eta from table
					call intptable(leta,lrhs,teta,trhs,NETA)

					eta(i,j) = teta

				end do
			end do

			return
		end

c-----------------------------------------------------------------------