CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     DRIVER FOR PS63 AT GALACTIC PROBLEM
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  "GALACTIC PROBLEM"
C  E. HAIRER, S.P. NORSET, G. WANNER,
C  SOLVING ORDINARY DIFFERENTIAL EQUATIONS, VOL. 1,
C  SPRINGER-VERLAG, 1991, SECOND EDITION, SECTION II.16
C----------------------------------------------
C LINK dr_Galactic.f RKPS63.f
C------------------------------------------
      PROGRAM GALACTIC
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(ND=6,NS=5,NDNS=ND*NS)
C----------------------------------------------
C      ND      : MAXIMAL DIMENSION OF THE PROBLEM
C      NS      : MAXIMAL NUMBER OF INTERNAL STAGES OF
C                THE EXPLICIT RUNGE-KUTTA METHOD
C----------------------------------------------
      PARAMETER(LYINT=2*NDNS,LYSOL=ND,LY=ND,LA=NS*(NS-1)/2,LZINT=2*NDNS)
      
      DIMENSION Y(LY),YSOL(LYSOL),YINT(LYINT),A(LA),B(NS),C(NS)
      CHARACTER*20 NOMFIC
C----------------------------------------------
C      Y       : ARRAY FOR THE NUMERICAL SOLUTION AT A GIVEN TIME
C      YSOL    : ARRAY FOR THE EXACT SOLUTION AT A GIVEN TIME
C      YINT    : ARRAY FOR WORKING SPACE
C      A       : ARRAY FOR THE MATRIX MA OF THE EXPLICIT RUNGE-KUTTA METHODS
C                FOR I>J, MA(I,J)=A((I-1)(I-2)/2+J)
C      B,C     : ARRAYS FOR THE VECTORS B AND C OF THE EXPLICIT RUNGE-KUTTA METHOD
C      NOMFIC  : NAME OF THE RESULT FILE
C----------------------------------------------
      EXTERNAL FEX,HAMEX,SOLOUT,SOLEX
      EXTERNAL INTEXP,RKPS63
C----------------------------------------------
C      FEX     : NAME OF SUBROUTINE COMPUTING THE VALUE OF F(T,Y)
C                OF THE FORM
C                    SUBROUTINE FEX(NDPB,T,Y,FY)    WHERE
C                NDPB  : DIMENSION OF Y AND F
C                T     : A GIVEN TIME
C                Y     : A GIVEN VECTOR
C                FY    : VALUE OF F(T,Y)
C      HAMEX   : NAME OF SUBROUTINE COMPUTING THE VALUE OF THE HAMILTONIAN H(T,Y)
C                OF THE FORM
C                    SUBROUTINE HAMEX(NDPB,T,Y,HAM)  WHERE
C                NDPB  : DIMENSION OF Y
C                T     : A GIVEN TIME
C                Y     : A GIVEN VECTOR
C                HAM   : VALUE OF H(T,Y)
C      SOLOUT  : NAME OF SUBROUTINE PROVIDING THE NUMERICAL SOLUTION DURING
C                INTEGRATION OF THE FORM
C                    SUBROUTINE SOLOUT(NDPB,T,Y,HAM,ERR,NFIC,NOINT)  WHERE
C                NDPB  : DIMENSION OF Y
C                T     : TIME OF THE OUTPUT
C                Y     : NUMERICAL SOLUTION AT TIME T
C                HAM   : VALUE OF THE HAMILTONIAN H(T,Y)
C                ERR   : EUCLIDEAN NORM OF THE VECTOR ERROR : Y_NUMERIC-Y_EXACT
C                        (IF Y_EXACT IS KNOWN)
C                NFIC  : NUMBER OF THE RESULT FILE
C                NOINT : NUMBER OF PERFORMED INTEGRATIONS AT TIME T
C      SOLEX   : NAME OF SUBROUTINE COMPUTING THE EXACT SOLUTION
C                (IF THERE EXISTS AN EXPLICIT FORM) OF THE FORM
C                    SUBROUTINE SOLEX(NDPB,T) WHERE
C                NDPB  : DIMENSION OF Y
C                T     : A GIVEN TIME
C----------------------------------------------
C      INTEXP  : NAME OF INTEGRATOR SUBROUTINE
C                OF THE FORM 
C                    SUBROUTINE INTEXP(NDPB,H,NBINT,T0,Y,HAM,FEX,HAMEX,
C                           SOLOUT,SOLEX,A,B,C,NFIC)  WHERE
C                NDPB   : DIMENSION OF Y AND F
C                H      : STEP SIZE
C                NBINT  : NUMBER OF INTEGRATIONS
C                T0     : INITIAL TIME
C                        REMARK : TEND=T0+NBINT*H
C                Y      : INPUT  : INITIAL VALUE OF Y
C                        OUTPUT : NUMERICAL SOLUTION AT TIME TEND
C                HAM    : INPUT  : INITIAL VALUE OF THE HAMILTONIAN
C                        OUTPUT : NUMERICAL HAMILTONIAN AT TIME TEND
C                FEX    : NAME OF SUBROUTINE COMPUTING F
C                HAMEX  : NAME OF SUBROUTINE COMPUTING H
C                SOLOUT : NAME OF SUBROUTINE PROVIDING THE NUMERICAL
C                         SOLUTION DURING INTEGRATION
C                SOLEX  : NAME OF SUBROUTINE COMPUTING EXACT SOLUTION
C                NFIC   : NUMBER OF THE RESULT FILE                
C      RKPS63  : NAME OF SUBROUTINE DEFINING THE EXPLICIT
C                RUNGE-KUTTA METHOD OF ORDER 3 AND
C                PSEUDO-SYMPLECTICNESS ORDER 6 WITH 5 STAGES,
C                CALLED PS63.
C                IT IS OF THE FORM
C                     SUBROUTINE RK-NAME-OF-THE-METHOD(A,B,C) WHERE
C                A   : ARRAY FOR THE MATRIX MA OF THE METHOD
C                      FOR I>J, MA(I,J)=A((I-1)(I-2)/2+J)
C                B,C : ARRAYS FOR THE VECTORS B AND C OF THE METHOD
C----------------------------------------------
	COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2
        COMMON/METHEXP/NSEXP
        COMMON/TRAVAIL/YINT
	COMMON/SORTIE/EXACT,YSOL
	COMMON/STAT/NFCN
C----------------------------------------------
C       /GALAXY/   : CONSTANTS FOR THE "GALACTIC PROBLEM"
C       /METHEXP/  : NSEXP  : DIMENSION OF THE VECTORS B AND C
C                             OF THE METHOD PS63
C	                    (NUMBER OF INTERNAL STAGES) 
C       /TRAVAIL/  : YINT   : ARRAY OF WORKING SPACE
C       /SORTIE/   : EXACT  : IF THERE EXISTS AN EXACT SOLUTION THEN
C                             EXACT<=1 ELSE EXACT>1
C                    YSOL   : VECTOR Y=(P(T),Q(T)) OF THE EXACT SOLUTION
C                             AT A GIVEN TIME T
C       /STAT/     : NFCN   : NUMBER OF F EVALUATIONS
C----------------------------------------------
C --- DEFINITION OF THE RESULT FILE
	WRITE(*,*)'NAME OF THE RESULT FILE?'
	READ(*,*) NOMFIC
	WRITE(*,*)'THE RESULT FILE IS : ',NOMFIC
	NFIC = 10
	OPEN(NFIC,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFIC)
	WRITE(NFIC,101)
 101	FORMAT('#',9X,'T',16X,' P',16X,' Q',16X,' HAM',16X,'ERR',10X,'NFCN')
C --- INITIALIZATION OF CONSTANTS SPECIFIC TO THE GALACTIC PROBLEM
	PAG    =  1.25D0
        PBG    =  1.0D0
        PCG    =  0.75D0
	PAG2   =  PAG*PAG
	PBG2   =  PBG*PBG
	PCG2   =  PCG*PCG
        AG     =  1.0D0
        CG     =  1.0D0
        OMEGAG =  0.25D0
C --- DIMENSION OF THE SYSTEM
        NDPB  =  6
C --- INITIAL VALUES
        T0     =  0.0D0
        Y(1)  =  0.0D0
        Y(2)  =  1.688837005904475591
        Y(3)  =  0.2D0
        Y(4)  =  2.5D0
        Y(5)  =  0.0D0
        Y(6)  =  0.0D0
C --- INITIAL HAMILTONIAN ("EXACT VALUE" OF THE HAMILTONIAN)
        HAM0  = 2.0D0
C --- THERE EXISTS NO EXPLICIT FORMULATION OF THE EXACT SOLUTION
C --- HENCE EXACT>1
	EXACT = 2.0D0
C --- INITIAL STEP SIZE
	WRITE(*,*)'STEP SIZE H? '
        READ(*,*) H
C --- NUMBER OF INTEGRATION NBINT OF THE FORM 
	WRITE(*,*)'NUMBER OF INTEGRATION? '
	READ(*,*) NBINT
C --- CALL OF THE SUBROUTINE RKPS63
      CALL RKPS63(A,B,C)
C --- CALL OF THE SUBROUTINE INTEXP
      CALL INTEXP(NDPB,H,NBINT,T0,Y,HAM0,FEX,HAMEX,SOLOUT,SOLEX,
     +            A,B,C,NFIC)
C --- CLOSE THE RESULT FILE 
      CLOSE(NFIC)
      STOP
      END

C----------------------------------------------
	SUBROUTINE FEX(NDPB,T,Y,FY)
C----------------------------------------------
C INPUT PARAMETERS
C       NDPB : DIMENSION OF THE PROBLEM (OF Y AND FY)
C       T    : A GIVEN TIME
C       Y    : A GIVEN VECTOR
C OUTPUT PARAMETERS
C       FY   : VALUE OF F(T,Y)
C----------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION Y(NDPB),FY(NDPB)
	COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2

	P1    =  Y(1)
        P2    =  Y(2)
        P3    =  Y(3)
	Q1    =  Y(4)
        Q2    =  Y(5)
        Q3    =  Y(6)
	QNR   =  CG + Q1*Q1/PAG2 + Q2*Q2/PBG2 + Q3*Q3/PCG2
        RIQNR =  1.0D0/QNR
        PNR   =  P1*P1 + P2*P2 + P3*P3
	FY(1) =  OMEGAG*P2 - 2.0D0*AG*Q1*RIQNR/PAG2
	FY(2) = -OMEGAG*P1 - 2.0D0*AG*Q2*RIQNR/PBG2
	FY(3) = -2.0D0*AG*Q3*RIQNR/PCG2
	FY(4) =  P1 + OMEGAG*Q2
	FY(5) =  P2 - OMEGAG*Q1
	FY(6) =  P3
	RETURN
        END

C----------------------------------------------
	SUBROUTINE HAMEX(NDPB,T,Y,HAM)
C----------------------------------------------
C INPUT PARAMETERS
C       NDPB : DIMENSION OF THE PROBLEM (OF Y AND FY)
C       T    : A GIVEN TIME
C       Y    : A GIVEN VECTOR
C OUTPUT PARAMETERS
C       HAM  : VALUE OF THE HAMILTONIAN H(T,Y)
C----------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION Y(NDPB)
	COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2

	P1    =  Y(1)
        P2    =  Y(2)
        P3    =  Y(3)
	Q1    =  Y(4)
        Q2    =  Y(5)
        Q3    =  Y(6)
	QNR   =  CG + Q1*Q1/PAG2 + Q2*Q2/PBG2 + Q3*Q3/PCG2
        RIQNR =  1.0D0/QNR
        PNR   =  P1*P1 + P2*P2 + P3*P3
	HAM   =  0.50D0*PNR + OMEGAG*(P1*Q2-P2*Q1) + AG*DLOG(QNR)
	RETURN
        END

C----------------------------------------------
	SUBROUTINE SOLEX(NDPB,T)
C----------------------------------------------
C INPUT PARAMETERS
C       NDPB : DIMENSION OF THE PROBLEM
C       T    : A GIVEN TIME
C OUTPUT PARAMETERS
C       YSOL : EXACT SOLUTION AT TIME T (IF THERE EXISTS AN
C              EXPLICIT FORM)
C----------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2
	COMMON/SORTIE/EXACT,YSOL(6)

	WRITE(*,*)'NO EXPLICIT FORMULATION FOR THE EXACT SOLUTION'
	RETURN
        END

C----------------------------------------------
 	SUBROUTINE SOLOUT(NDPB,T,Y,HAM,ERROR,NFIC,NOINT)
C----------------------------------------------
C INPUT PARAMETERS
C     NDPB  : DIMENSION OF Y
C     T     : TIME OF THE OUTPUT
C     Y     : NUMERICAL SOLUTION AT TIME T
C     HAM   : VALUE OF THE HAMILTONIAN H(T,Y)
C     ERROR : EUCLIDEAN NORM OF THE VECTOR ERROR : Y_NUMERIC-Y_EXACT
C                        (IF Y_EXACT IS KNOWN)
C    NFIC  : NUMBER OF THE RESULT FILE
C    NOINT : NUMBER OF PERFORMED INTEGRATIONS AT TIME T
C PURPOSE : WRITE ON THE RESULT FILE NUMBER "NFIC" THE TIME T, TWO
C           CHOSEN COMPONENTS OF THE NUMERICAL SOLUTION AT TIME T,
C           THE NUMERICAL HAMILTONIAN, THE ERROR AND THE NUMBER OF
C           F EVALUATIONS
C----------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        COMMON/STAT/NFCN
	DIMENSION Y(NDPB)

C --- FOR AN OUTPUT AT ALL THE 100 INTEGRATIONS :
        NR = MOD(NOINT,100)
        IF((NOINT.EQ.0).OR.(NR.EQ.0))THEN
  	WRITE(NFIC,110)T,Y(1),Y(2),HAM,ERROR,NFCN
        ENDIF
 110	FORMAT(1X,5E18.10,I10)
	RETURN
        END
