!---------------------------------------------------------------------- 
!	sphere.f 
!	Contains routines to construct an initial surface  
!	Modified 5-21-93	MFH 
!		- Switched to Spherical coordinates 
!		- S(theta,phi,3) -> 1 -> radial, 2 -> theta, 3-> phi 
!		- Indices for theta-phi mesh on surface 
!	Modified 6-30-93	FH 
!		- Switched to non-staggered grid 
!	Modified Sat Sep 18 18:14:27 CDT 1993 MFH 
!		- Center of coordinate system center(3) 
!---------------------------------------------------------------------- 
 
	subroutine GenSur(S,Sc,Sf,NDIM,NST,NSP) 
		implicit none 
		include 'parameters.h' 
 
		integer NDIM 
		integer NST 
		integer NSP 
		integer np,nt 
		 
		real S(NST,NSP,NDIM) 
		real Sc(NST,NSP,NDIM) 
		real Sf(NST,NSP) 
 
 
 
!	Local variables 
 
		integer i,j,k 
		 
		real x,y,z 
		real t,p 
		real temp 
 
!		write(6,*)'Radius ',radius 
!		write(6,*)'Number ',NSS*(NSS+1) 
!		write(0,*)x0,y0,z0,tau 
		write(0,*)'CENTER',center(1),center(2),center(3) 
		write(2,*)'CENTER',center(1),center(2),center(3) 
 
!	Define PI 
		Pi = 4.0e0*atan(1.0e0) 
 
!	Do loop to calculate points on a sphere 
 
!	Store values as r,theta,phi 
		do i=1,NST 
			t = float(i-1)*dtheta 
			do j=1,NSP 
				 p = float(j-1)*dphi 
 
				 x = radius*cos(p)*sin(t) + center(1) 
				 y = radius*sin(p)*sin(t) + center(2) 
				 z = radius*cos(t)        + center(3) 
 
!				 Sf(i,j)  = 0.0e0 
!				 Sf(i,j)  = sqrt(x**2 + y**2 + z**2) 
!				 Sf(i,j)  = radius 
!				 S(i,j,1) = radius 
 
				 Sf(i,j)  = sqrt((x-center(1))**2 + (y-center(2))**2 + 
     #								 (z-center(3))**2) 
!			Spherical coordinates 
!				 S(i,j,1)  = sqrt((x-center(1))**2 + (y-center(2))**2 + 
!     #								 (z-center(3))**2) 
				 S(i,j,1) = Sf(i,j) 
				 S(i,j,2) = t 
				 S(i,j,3) = p 
 
!			Cartesian Coordinates 
				 Sc(i,j,1) = x 
				 Sc(i,j,2) = y 
				 Sc(i,j,3) = z 
 
				 
!				 write(66,*)S(i,j,1),S(i,j,2),S(i,j,3) 
!				 write(67,*)Sc(i,j,1),Sc(i,j,2),Sc(i,j,3) 
!				 write(68,*)x,y,z 
!				 write(3,*)i,j,Sf(i,j) 
			enddo 
!			write(3,*) 
!			write(68,*) 
		enddo 
		return 
		end 
!---------------------------------------------------------------------- 
!	Modified routine that does not require Sf	 
! 
	subroutine InitSur(S,Sc,NST,NSP) 
		implicit none 
		include 'parameters.h' 
 
		integer NST 
		integer NSP 
		integer np,nt 
		 
		real S(NST,NSP,3) 
		real Sc(NST,NSP,3) 
 
 
 
!	Local variables 
 
		integer i,j,k 
		 
		real x,y,z 
		real t,p 
		real temp 
 
		write(0,*)'CENTER',center(1),center(2),center(3) 
		write(2,*)'CENTER',center(1),center(2),center(3) 
 
!	Define PI 
		Pi = 4.0e0*atan(1.0e0) 
 
!	Do loop to calculate points on a sphere 
 
!	Store values as r,theta,phi 
		do i=1,NST 
			t = float(i-1)*dtheta 
			do j=1,NSP 
				 p = float(j-1)*dphi 
 
				 x = radius*cos(p)*sin(t) + center(1) 
				 y = radius*sin(p)*sin(t) + center(2) 
				 z = radius*cos(t)        + center(3) 
 
				 S(i,j,1) = sqrt((x-center(1))**2 + (y-center(2))**2 + 
     #                         (z-center(3))**2) 
				 S(i,j,2) = t 
				 S(i,j,3) = p 
 
!			Cartesian Coordinates 
				 Sc(i,j,1) = x 
				 Sc(i,j,2) = y 
				 Sc(i,j,3) = z 
 
				 
			enddo 
		enddo 
		return 
		end 
!---------------------------------------------------------------------- 
