MODULE values

  IMPLICIT NONE
  
  REAL*8, PARAMETER :: omegan   = 0.21D0
  REAL*8, PARAMETER :: radk     = 12.D0
  REAL*8, PARAMETER :: a        = 0.998D0
  REAL*8, PARAMETER :: nd_c     = 1.D11
  REAL*8, PARAMETER :: gbeta    = 1.D-05
  REAL*8, PARAMETER :: inkappa  = 4.301D12 
  REAL*8, PARAMETER :: Mass     = 1.D0
  REAL*8, PARAMETER :: piby2    = 2.D0*DATAN(1.D0)
  
  INTEGER, PARAMETER :: NEQ = 1
  INTEGER, PARAMETER :: Nmax = 200000
  
END MODULE values


MODULE direction
   
  IMPLICIT NONE
  
  INTEGER :: dirn
  INTEGER :: num_eq
  REAL*8 :: gbeta13, gbeta43
  REAL*8 :: r_val
  
  REAL*8, POINTER :: r_eq(:)
  REAL*8, POINTER :: nd_eq(:)
  INTEGER, POINTER :: num_nd(:)
  
END MODULE direction


MODULE Calculations

  USE values
  USE direction
  
  CONTAINS
  
  
  FUNCTION sigma(r,theta)
  
    IMPLICIT NONE
	
	REAL*8 :: sigma, r, theta
	
	sigma = r*r + (a*DCOS(theta))**2
	
  END FUNCTION sigma
  
  
  FUNCTION delta(r)
  
    IMPLICIT NONE
	
	REAL*8 :: delta, r
	
	delta = r*r - 2.D0*Mass*r + a*a
	
  END FUNCTION delta
  
  
  FUNCTION omega(x)
  ! Compute the angular velocity law \Omega for the torus
    IMPLICIT NONE
	
	REAL*8 :: omega, x
! There should be an M-dependence here, add it later!	
	omega = ((radk/x)**omegan)/(a+x*DSQRT(x))
	
  END FUNCTION omega
  
  
  FUNCTION td2(r,theta)
  ! Calculate tdot^2
    IMPLICIT NONE
	
	REAL*8 :: td2, r, theta
	REAL*8 :: siny, siny2, om, tmr, sig, del, insig
	
	siny  = DSIN(theta)
	siny2 = siny*siny
	om    = omega(r*siny)
	tmr   = 2.D0*Mass*r
	sig   = sigma(r,theta)
	del   = delta(r)
	insig = 1.D0/sig
	
	td2 = insig*(sig-tmr+2.D0*tmr*a*siny2*om-(sig*del+tmr*(r*r+a*a))*siny2*om*om)
	td2 = 1.D0/td2
	
  END FUNCTION td2
  
  
  FUNCTION r_density(N1,X,Y)
  ! ODE of the density of the torus in the equatorial plane
    IMPLICIT NONE
	
	REAL*8  :: r_density
	INTEGER :: N1
	REAL*8  :: X
	REAL*8, DIMENSION(NEQ) :: Y
	REAL*8  :: inX, om, temp
	
	inX  = 1.D0/X
	om   = omega(X)
	temp = (1.D0-a*om)*inX
	
	r_density = 3.D0*Y(1)*(X*om*om-Mass*temp*temp)*&
	(1.D0+inkappa*gbeta43*Y(1)**(-1.D0/3.D0))*td2(X,piby2)
	
  END FUNCTION r_density
  
  
  SUBROUTINE fireequator(yn,counter)
  ! Integrate along the equatorial plane
    IMPLICIT NONE
	
	REAL*8, DIMENSION(NEQ) :: Y0, YNN, Yf
	REAL*8, DIMENSION(NEQ,12) :: WORK
	REAL*8, DIMENSION(Nmax) :: XT, F1
	REAL*8, DIMENSION(2) :: yn
	REAL*8 :: X0, XN, DX, EPS
	INTEGER :: ierr, i, counter
	
	EPS = 1.D-09
	XT(1) = yn(1)
	counter = 0
	DX = 1.D-04
	F1(1) = yn(2)

	
	DO i = 1, Nmax-1
	
	  IF (F1(i) < 1.D0) THEN
	    DX = 1.D-06
	  ELSE IF (F1(i) < 1.D02) THEN
	    DX = 1.D-05
	  END IF
	  
	  XT(i+1) = XT(i) + DBLE(dirn)*DX
	  X0 = XT(i)
	  Y0(1) = F1(i)
	  XN = XT(i+1)
	  
	  CALL RUNGE(NEQ,r_density,X0,XN,Y0,EPS,YNN,ierr,WORK)
	  
	  F1(i+1) = YNN(1)
	  
	  IF (.NOT. (F1(i+1)<=0.D0 .OR. F1(i+1)>=0.D0)) THEN
	    EXIT 
	  END IF
	  
	  counter = counter + 1
!      Write Log10(rho)	  
!	   WRITE(10,'(2ES25.15)') xt(i+1), DLOG10(F1(i+1))
	  WRITE(10,'(2ES25.15)') xt(i+1), F1(i+1)
	  
	END DO
	
  END SUBROUTINE fireequator
  
  
  SUBROUTINE calc_equator(n1,n2)
  
    IMPLICIT NONE
	
	INTEGER :: n1, n2, counter
	REAL*8, DIMENSION(2) :: yn
	
	OPEN(UNIT=10, FILE='equator.dat', STATUS='OLD', ACTION='WRITE')
	
	yn(1) = radk   ! radius, r
	yn(2) = nd_c   ! density, rho
	dirn  = -1
!	 WRITE(10,'(2ES25.15)') yn(1), DLOG10(yn(2))
    WRITE(10,'(2ES25.15)') yn(1), yn(2)
	CALL fireequator(yn,counter)
	
	n1 = counter + 1
	
	yn(1) = radk
	yn(2) = nd_c
	dirn  = 1
	CALL fireequator(yn,n2)
	
	CLOSE(10)
	
  END SUBROUTINE calc_equator
  
  
  SUBROUTINE read_equator (num_eqat1, num_eqat2)

    IMPLICIT NONE

    INTEGER :: num_eqat1, num_eqat2
    INTEGER :: i

    OPEN(UNIT=10, FILE='equator.dat',STATUS='OLD',ACTION='READ')

    ! Backward direction first
    DO i = 1, num_eqat1
      READ (10,'(2ES25.15)') r_eq(num_eqat1+1-i), nd_eq(num_eqat1+1-i)
    END DO

    ! Now the forward section

    DO i = 1, num_eqat2
      READ (10,'(2ES25.15)') r_eq(num_eqat1+i), nd_eq(num_eqat1+i)
    END DO

    CLOSE(UNIT=10)

  END SUBROUTINE read_equator
  
  
  SUBROUTINE ordered_equatorial()
  ! Order the (r,rho) values from low to high r  
    IMPLICIT NONE
	
	INTEGER :: i
	
	OPEN (UNIT=11,FILE='order_equator.dat',STATUS='OLD',ACTION='WRITE')
	
	DO i = 1, num_eq
	  WRITE(11,'(2ES25.15)') r_eq(i), nd_eq(i)
	END DO
	
	CLOSE(UNIT=11)
	
  END SUBROUTINE ordered_equatorial
  
  
  
  FUNCTION theta_density(N2,X,Y)
  ! ODE of the density of the torus as a function of theta for fixed r
    IMPLICIT NONE
	
	REAL*8  :: theta_density
	INTEGER :: N2
	REAL*8  :: X
	REAL*8, DIMENSION(NEQ) :: Y
	
	REAL*8  :: r, theta, siny, cosy, om
	REAL*8  :: sig, del, temp 

    r    = r_val
	theta= DASIN(X)
	siny = X
	cosy = DSQRT(1.D0-siny*siny)
	om   = omega(r*siny)
	
	sig  = sigma(r,theta)
	del  = delta(r)
	temp = (a-(r*r+a*a)*om)/sig
	
	theta_density = 3.D0*Y(1)*siny*cosy*(2.D0*Mass*r*temp*temp+del*om*om)*&
	                (1.D0+inkappa*gbeta43*Y(1)**(-1.D0/3.D0))*td2(r,theta)
					
  END FUNCTION theta_density
  
  
  SUBROUTINE firerho(yn,ncount)
  ! Integrate along constant-r surfaces for varying theta until rho=0
    IMPLICIT NONE
	
	REAL*8, DIMENSION(NEQ) :: Y0, YNN, Yf
	REAL*8, DIMENSION(NEQ,12) :: WORK
	REAL*8, DIMENSION(Nmax) :: XT, F1
	REAL*8, DIMENSION(2) :: yn
	REAL*8 :: X0, XN, DX, EPS
	INTEGER :: ierr, i, ncount
	
	EPS = 1.D-09
	XT(1) = yn(1)  ! SIN(theta)
	ncount = 0
	DX = 1.D-04
	F1(1) = yn(2)  ! rho

	WRITE(5,'(2ES25.15)') yn(1), yn(2)
	
	DO i = 1, Nmax-1
	
	  IF (F1(i) < 1.D0) THEN
	    DX = 1.D-06
	  ELSE IF (F1(i) < 1.D02) THEN
	    DX = 1.D-05
	  END IF
	  
	  XT(i+1) = XT(i) - DX
	  X0 = XT(i)
	  Y0(1) = F1(i)
	  XN = XT(i+1)
	  
	  CALL RUNGE(NEQ,theta_density,X0,XN,Y0,EPS,YNN,ierr,WORK)
	  
	  F1(i+1) = YNN(1)
	  
	  IF (.NOT. (F1(i+1)<=0.D0 .OR. F1(i+1)>=0.D0)) THEN
	    EXIT 
	  END IF
	  
	  IF (F1(i+1) <= 0.D0) THEN
	    EXIT
	  END IF
	  
	 ncount = ncount + 1
!      Write Log10(rho)	  
!	   WRITE(5,'(2ES25.15)') xt(i+1), DLOG10(F1(i+1))
	  WRITE(5,'(2ES25.15)') xt(i+1), F1(i+1)
	  
	END DO
	
  END SUBROUTINE firerho
  
  
  SUBROUTINE calc_density(maxn)
  
    IMPLICIT NONE
	
	INTEGER :: maxn, i
	REAL*8, DIMENSION(2) :: yn
	
	OPEN(UNIT=5,FILE='torus.dat',STATUS='OLD',ACTION='WRITE')
	
	maxn = 0
	
	DO i = 1, num_eq
	  yn(1) = 1.D0
	  yn(2) = nd_eq(i)
	  r_val = r_eq(i)
	  WRITE(5,'(1F25.15)') r_val
	  Call firerho(yn,num_nd(i))
	  
	  IF (num_nd(i) > maxn) THEN
	    maxn = num_nd(i)
	  END IF
	  
	END DO
	
	CLOSE(UNIT=5)
	
  END SUBROUTINE calc_density
   
   
   
END MODULE Calculations


     SUBROUTINE RUNGE(NEQ,FUNC,X0,XN,Y0,EPS,YN,IERR,WORK) 
! ------------------------------------------------------------          
! SOLVE SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS                       
!                                                                       
!        DY(  1)/DX = FUNC(  1,X,Y)                                     
!        DY(  2)/DX = FUNC(  2,X,Y)                                     
!           ...............                                             
!        DY(NEQ)/DX = FUNC(NEQ,X,Y)                                     
!                                                                       
!       BY 5-TH ORDER RUNGE KUTTA METHOD BASED ON                       
!          FEHLBERG'S FORMULA IN DOUBLE                                 
!              IN DOUBLE PRECISION                                      
!                                                                       
!   THIS ROUTINE CALLS A 5-TH ORDER RUNGE-KUTTA TYPE FORMULA,           
!   BUT ERROR ESTIMATION IS CARRIED BY COMPARING THE SOLUTION           
!   WITH A 4-TH ORDER FORMULA.                                          
!                                                                       
!   ==== INPUT ====                                                     
!                                                                       
!   NEQ ... ORDER OF DIFFERENTIAL EQUATION                              
!   FUNC ... NAME OF FUNCTION SUBPROGRAM FOR RIGHT HAND SIDE            
!            OF THE EQUATIONS                                           
!     EXAMPLE : FUNC=EX2, NEQ=2, X=T, Y(1)=Y, Y(2)=-X                   
!                                                                       
!         DX/DT =  Y                                                    
!         DY/DT = -X                                                    
!                                                                       
!   X0 ... INITIAL VALUE OF INDEPENDENT VARIABLE X                      
!   XN ... VALUE OF INDEPENDENT VARIABLE X AT WHICH                     
!          APPROXIMATE SOLUTION IS REQUIRED                             
!   Y0(I),I=1,NEQ ... INITIAL VALUE AT X0                               
!       Y0(I) WILL BE CHANGED IN THIS SUBROUTINE                        
!   EPS ... ABSOLUTE ERROR TOLERANCE FOR YN                             
!                                                                       
!   ==== OUTPUT ====                                                    
!                                                                       
!   YN(I),I=1,NEQ ... APPROXIMATE SOLUTION AT XN                        
!   IERR ... IF IERR=0 THEN CONVERGED                                   
!            IF IERR=1 THEN NOT CONVERGED AND COMPUTION TERMINATED      
! -----------------------------------------------------------------     
      IMPLICIT REAL*8(A-H,O-Z) 
      DIMENSION Y0(NEQ),YN(NEQ),WORK(NEQ,12) 
!                                                                       
      DATA ITEMAX / 100 / 
      DATA EPSMAC / 1.D-16 / 
!                                                                       
      EXTERNAL FUNC 
!                                                                       
      IF(EPS >= EPSMAC) THEN 
        EPSV = EPS 
       ELSE 
        EPSV = EPSMAC 
      END IF 
!                                                                       
      EB = ABS( EPSV / (XN - X0) ) 
!                                                                       
!     ---- INITIAL STEP ----                                            
!                                                                       
      IERR = 0 
!                                                                       
      X = X0 
      H = XN - X0 
!                                                                       
      DO 10 ITER = 1, ITEMAX 
!                                                                       
!       ---- GENERAL STEP ----                                          
!                                                                       
        CALL FEHLBG (NEQ, FUNC, X, H, Y0, YN, ERRMAX,                   &
     &    WORK(1,1), WORK(1,2),  WORK(1,3),  WORK(1,4),                 &
     &    WORK(1,5), WORK(1,6),  WORK(1,7),  WORK(1,8),                 &
     &    WORK(1,9), WORK(1,10), WORK(1,11), WORK(1,12))                
!                                                                       
        IF (ERRMAX <= ABS(EB*H)) THEN 
          X = X + H 
          H = XN - X 
!                                                                       
!         ---- CONVERGENCE TEST ----                                    
!                                                                       
          IF (ABS(H) <= EPSMAC) RETURN 
!                                                                       
          DO 520 I = 1, NEQ 
            Y0(I) = YN(I) 
  520     CONTINUE 
!                                                                       
         ELSE 
          IF (ERRMAX > EPSV) THEN 
            H = 0.1*H 
           ELSE 
            H = 0.9*H*SQRT(SQRT(EB*ABS(H)/ERRMAX)) 
            IF (ABS(H) < EPSMAC) GO TO 900 
          END IF 
        END IF 
!                                                                       
   10 END DO 
!                                                                       
!     ---- NOT CONVERGED ----                                           
!                                                                       
  900 CONTINUE 
      IERR = 1 
!      WRITE(*,2001) X 
! 2001 FORMAT(' (RUNGE) TROUBLE AT X = ',1PD13.6,                        &
!     & '  RETURN WITHOUT FURTHER COMPUTATION')                          
      RETURN 
                                                                      
 END SUBROUTINE

                                          
                                                                     
 SUBROUTINE FEHLBG (NEQ, FUNC, X, H, Y, YN, ERRMAX,                &
     &         AK1, AK2, AK3, AK4, AK5, AK6,                            &
     &         W2,  W3,  W4,  W5,  W6,  ERR)                            
! -----------------------------------------------------------           
!        5-TH ORDER RUNGE-KUTTA FORMULA                                 
!               BY FEHLBERG                                             
!                                                                       
!     Y(I) ... INITIAL VALUE AT X                                       
!     YN(I) ... SOLUTION AT X + H                                       
!     ERRMAX ... ESTIAMTE OF MAX ERROR OF SOLUTION                      
! -----------------------------------------------------------           
      IMPLICIT REAL*8(A-H,O-Z) 
      DIMENSION Y(NEQ),YN(NEQ),ERR(NEQ) 
      DIMENSION AK1(NEQ),AK2(NEQ),AK3(NEQ),                             &
     &          AK4(NEQ),AK5(NEQ),AK6(NEQ)                              
      DIMENSION W2(NEQ),W3(NEQ),W4(NEQ),W5(NEQ),W6(NEQ) 
!                                                                       
      PARAMETER (ONE = 1, THR = 3, TWL = 12) 
!                                                                       
      PARAMETER (AL2 = ONE / 4,  AL3 = THR / 8,                         &
     &           AL4 = TWL / 13, AL5 = ONE,                             &
     &           AL6 = ONE / 2)                                         
!                                                                       
      PARAMETER (BD2 = 4, BD3 = 32, BD4 = 2197, BD5 = 4104,             &
     &           BD6 = 20520)                                           
!                                                                       
      PARAMETER (B21 =      1/BD2,                                      &
     &           B31 =      3/BD3, B32 =      9/BD3,                    &
     &           B41 =   1932/BD4, B42 =  -7200/BD4,                    &
     &           B43 =   7296/BD4,                                      &
     &           B51 =   8341/BD5, B52 = -32832/BD5,                    &
     &           B53 =  29440/BD5, B54 =   -845/BD5,                    &
     &           B61 =  -6080/BD6, B62 =  41040/BD6,                    &
     &           B63 = -28352/BD6, B64 =   9295/BD6,                    &
     &           B65 =  -5643/BD6)                                      
!                                                                       
      PARAMETER (GAD = 7618050) 
      PARAMETER (GA1 =  902880/GAD, GA3 =  3953664/GAD,                 &
     &           GA4 = 3855735/GAD, GA5 = -1371249/GAD,                 &
     &           GA6 =  277020/GAD)                                     
!                                                                       
      PARAMETER (DAD = 752400) 
      PARAMETER (DA1 =   2090/DAD, DA3 = -22528/DAD,                    &
     &           DA4 = -21970/DAD, DA5 =  15048/DAD,                    &
     &           DA6 =  27360/DAD)                                      
!                                                                       
      IF (NEQ > 10) RETURN 
!                                                                       
      DO I = 1, NEQ 
        AK1(I) = H*FUNC(I, X, Y) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        W2(I) = Y(I)+B21*AK1(I) 
      END DO 
!
      DO I = 1, NEQ 
        AK2(I) = H*FUNC(I, X+AL2*H, W2) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        W3(I) = Y(I)+B31*AK1(I)+B32*AK2(I) 
      END DO 
!
      DO I = 1, NEQ 
        AK3(I) = H*FUNC(I, X+AL3*H, W3) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        W4(I) = Y(I)+B41*AK1(I)+B42*AK2(I)+B43*AK3(I) 
      END DO 
!      
      DO I = 1, NEQ 
        AK4(I) = H*FUNC(I, X+AL4*H, W4) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        W5(I) = Y(I)+B51*AK1(I)+B52*AK2(I)+B53*AK3(I)+B54*AK4(I) 
      END DO
! 
      DO I = 1, NEQ 
        AK5(I) = H*FUNC(I, X+AL5*H, W5) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        W6(I) = Y(I)+B61*AK1(I)+B62*AK2(I)+B63*AK3(I)+B64*AK4(I)        &
     &              +B65*AK5(I)                                         
      END DO 

      DO I = 1, NEQ 
        AK6(I) = H*FUNC(I, X+AL6*H, W6) 
      END DO 
!                                                                       
      DO I = 1, NEQ 
        YN(I) = Y(I)+GA1*AK1(I)+GA3*AK3(I)+GA4*AK4(I)+GA5*AK5(I)        &
     &              +GA6*AK6(I)                                         
      END DO 
!                                                                       
      DO I = 1, NEQ 
        ERR(I) = DA1*AK1(I)+DA3*AK3(I)                                  &
     &          +DA4*AK4(I)+DA5*AK5(I)+DA6*AK6(I)                       
      END DO 
!                                                                       
      ERRMAX = 0.D0 
      DO I = 1, NEQ 
        IF( ABS(ERR(I)) .GT. ERRMAX) ERRMAX = ABS(ERR(I)) 
      END DO 
!                                                                       
      RETURN 

 END SUBROUTINE       


PROGRAM MAIN
  USE values
  USE direction
  USE Calculations
  IMPLICIT NONE
  
  INTEGER :: n1, n2, maxn
  
  gbeta13 = gbeta**(1.D0/3.D0)
  gbeta43 = gbeta**(4.D0/3.D0)
  
  CALL calc_equator(n1,n2)
  
  num_eq = n1 + n2

  ALLOCATE(r_eq(num_eq))
  ALLOCATE(nd_eq(num_eq))
  
  CALL read_equator(n1,n2)
  
!  CALL ordered_equatorial

  ALLOCATE(num_nd(num_eq))
! Remember to deallocate memory when no-longer needed  
  WRITE(*,*) 'The equatorial range is'
  WRITE(*,'(2F25.15)') r_eq(1), r_eq(num_eq)
  
  CALL calc_density(maxn)
  
END PROGRAM MAIN
