PROGRAM SPAM25
C
C       COPYRIGHT 1987 STARCHILD SOFTWARE
C
C      GALAXY INTERACTIONS
C      AUTHOR: JOHN F. WALLIN
C     THIS PROGRAM READS IN A FILE WITH THE INITIAL
C     TEST POINT POSITIONS AND VELOCITIES ALONG WITH THE
C     GALAXY POTENTIAL CHARACTERISTICS AND EVOLVES THEM.
C     THE DATA FILE REQUIRED IS OF THE FORM-
C      T0  = INITIAL TIME
C       N  = NUMBER OF TEST POINTS + 1
C     X(I),VX(I),Y(I),VY(I),Z(I),VZ(I)
C      POSITIONS AND VELOCITIES OF THE N-1 TEST PARTICLES
C      AND POSITIONS AND VELOCITIES OF THE INTERSECTING
C      GALAXY.  THE MAIN GALAXY IS AT (0,0,0) WITH V=0.
C      CP(1-6)
C      CM(1-6)
C      THE CP AND CM PARAMETERS ARE USED TO CHARACTERIZE
C      THE POTENTIALS OF THE PERTURBING AND MAIN  GALAXIES.
C      THEY ARE USED IN THE ROUTINE DIFFEQ TO FIND THE
C      GRAVITATIONAL ACC.
C      THE OUTPUT FILES ARE STORED IN THE SAME FORMAT AS
C      THE INPUT FILE.
C      FOR SPECIFIC FORMATS USED, SEE OUTPUT AND INPUT SUBROUTINES.
C
C
       DIMENSION X0(10001,7),XE(10001,7),CP(6),CM(6),FRICTN(6)
       DIMENSION XPOS(10001),YPOS(10001)
       REAL*8 X0,XE,CP,CM,TF,DT,P0,H,ALPHA,BETA
       REAL*8 INCL,OMEGA,RMIN,VCAR
       REAL*8 T1(10),T2(2)
        REAL*8 RSCALE,XSIZE,YSIZE
        INTEGER NF,III
        INTEGER I,J,N
        INTEGER XPOS,YPOS,GRID(200,200),NN
        CHARACTER*10 NAME
	CHARACTER*60 OUTFILE10, INFILE
	INTEGER*2 FORM1
        INTEGER N1,N2,NPTS,PT(60),NUM(60)
C
        COMMON /POS/ X0,XE,CM,CP,N
        COMMON /PROJ/ XPOS,YPOS,RSCALE,XSIZE,YSIZE
        COMMON /INP/ N1,N2,H,ALPHA,BETA
        COMMON /FISH/ GRID
C
	PRINT*,'INPUT FILE NAME?'
	READ(*,'(A)') INFILE
	PRINT*,'POSITION/VELOCITY OUTPUT FILE NAME?'
	READ(*,'(A)') OUTFILE10
	PRINT*,'OUTPUT FILE FORM:'
	PRINT*,'0. ASCII'
	PRINT*,'1. COMPRESSED BINARY'
	READ*,FORM1
C
	PRINT*,'PRINT FINAL TIME?'
	READ*,TF
	PRINT*,'OUTPUT TIME STEP?'
	READ*,DT
        OPEN(11,FILE=INFILE,STATUS='OLD')
	CALL INPUT
	PRINT*,'DATA READ'
C
	IF (FORM1.EQ.0) THEN
	OPEN(10,FILE=OUTFILE10, STATUS='NEW')
	ELSE
	OPEN(10,FILE=OUTFILE10,FORM='UNFORMATTED',STATUS='NEW')
	ENDIF

C      END OF PROGRAM ?
40     IF (TF.LT.X0(1,1)) GOTO 10

C       PRINT OUTPUT ??
       P0=X0(1,1)
       P0=ABS(P0) + ABS(H/2.)
       P0 = DMOD(P0,DT)

        IF (P0.LT.H) THEN
	  IF (FORM1.EQ.0) THEN 
            CALL OUTPT2	
	  ELSE
            CALL OUT2B
	  ENDIF 
	   dist = sqrt(x0(n,2)**2 + x0(n,4)**2 + x0(n,6)**2)
           PRINT*,'TIME = ',X0(1,1),' dist =',dist      
         ENDIF

C      MOVE STUFF
        CALL RK4(H)

C      REDEFINE X0
         DO 30 I=1,N
         DO 30 J=1,7
30       X0(I,J)=XE(I,J)
C
C      REPEAT PROCEDURE
         GOTO 40
C
C
10      CONTINUE
	  IF (FORM1.EQ.0) THEN 
            CALL OUTPT2	
	  ELSE
            CALL OUT2B
	  ENDIF 
        CLOSE(10)
        CLOSE(11)
        PRINT*,'END'
        STOP
        END

C------------------------------------------------------------------
C
       SUBROUTINE OUTPUT
*
       DIMENSION X0(10001,7),CP(6),CM(6)
       DIMENSION XE(10001,7),XPOS(10001),YPOS(10001)
       REAL*8 CP,CM,XE,X0,H,ALPHA,SALPHA,CALPHA,PI,XT,YT
       REAL*8 SBETA,CBETA,BETA
       REAL*8 RSCALE,XSIZE,YSIZE
       INTEGER N,I
       INTEGER XPOS,YPOS
       INTEGER N1,N2
       COMMON /POS/ X0,XE,CM,CP,N
       COMMON /PROJ/ XPOS,YPOS,RSCALE,XSIZE,YSIZE
       COMMON /INP/ N1,N2,H,ALPHA,BETA
C
        PI=3.1415926
C
        SALPHA=SIN(ALPHA*PI/180.)
        CALPHA=COS(ALPHA*PI/180.)
        SBETA=SIN(BETA*PI/180.)
        CBETA=COS(BETA*PI/180.)
C
C
        DO 10 I=1,N
        XT=X0(I,2)*CBETA-X0(I,4)*SBETA
        YT=X0(I,2)*SBETA+X0(I,4)*CBETA
        XPOS(I)=INT(XT*RSCALE+XSIZE/2.)
        YPOS(I)=INT((YT*CALPHA-X0(I,6)*SALPHA)*RSCALE+YSIZE/2.)
10      CONTINUE
C
C
        RETURN
        END
C
C*
C-----------------------------------------------------------------

        SUBROUTINE INPUT
        DIMENSION X0(10001,7),XE(10001,7),CP(6),CM(6)
        REAL*8 X0,XE,CP,CM,ALPHA,BETA,H
        INTEGER N,I,N1,N2
        INTEGER NNX,NN1,NN2
        INTEGER NF
        CHARACTER TMP*10
       COMMON /POS/ X0,XE,CM,CP,N
       COMMON /INP/ N1,N2,H,ALPHA,BETA
C
*
*
*       LOCAL GRID SIZE
        NN = 100
*
*       OUTPUT INTERVAL AND FINISH TIME
*        DT = 1.
*        TF = 1.0
*
       READ(11,2) N, N1, N2
        PRINT*,'N = ',N
       READ(11,4) NN1,NN2
       READ(11,4) NN1,NN2
       READ(11,1) ALPHA,BETA,H,X0(1,1)
C
       DO 10 I=1,N
       READ(11,3) X0(I,2),X0(I,4),X0(I,6),X0(I,3),X0(I,5),X0(I,7)
10     CONTINUE
       READ(11,3) CP(1),CP(2),CP(3),CP(4),CP(5),CP(6)
       READ(11,3) CM(1),CM(2),CM(3),CM(4),CM(5),CM(6)
       READ(11,5) TMP
C
C
1     FORMAT(4E16.8)
2     FORMAT(3I5)
3     FORMAT(6E16.8)
4     FORMAT(2I5)
5      FORMAT(A10)
*
*
        RETURN
        END
C
C
C-----------------------------------------------------------------
        SUBROUTINE OUTPT2
       DIMENSION X0(10001,7),XE(10001,7),CP(6),CM(6)
        REAL*8 X0,XE,CP,CM,ALPHA,BETA,H
        INTEGER N,I,N1,N2
        CHARACTER NAME*10
       COMMON /POS/ X0,XE,CM,CP,N
       COMMON /INP/ N1,N2,H,ALPHA,BETA
C
       WRITE(10,2) N,N1,N2
       WRITE(10,4) 0,0
       WRITE(10,4) 0,0
       WRITE(10,1) ALPHA,BETA,H,X0(1,1)
C
       DO 10 I=1,N
       WRITE(10,3) X0(I,2),X0(I,4),X0(I,6),X0(I,3),X0(I,5),X0(I,7)
10     CONTINUE
       WRITE(10,3) CP(1),CP(2),CP(3),CP(4),CP(5),CP(6)
       WRITE(10,3) CM(1),CM(2),CM(3),CM(4),CM(5),CM(6)
       WRITE(10,5)
C
C
 1     FORMAT(4E16.8)
 2     FORMAT(3I5)
 3     FORMAT(6E16.8)
 4     FORMAT(2I5)
 5     FORMAT(1X,'END')
        RETURN
        END
C
C

C-----------------------------------------------------------------
        SUBROUTINE OUT2B
       DIMENSION X0(10001,7),XE(10001,7),CP(6),CM(6)
        REAL*8 X0,XE,CP,CM,ALPHA,BETA,H
        INTEGER N,I,N1,N2
        CHARACTER NAME*10
	INTEGER*2 NN2,N21,N22
	INTEGER*2 ALPH,BET,H2,T2
	INTEGER*2 X(6),CP2(6),CM2(6)
 
	REAL XZZ
	
       COMMON /POS/ X0,XE,CM,CP,N
       COMMON /INP/ N1,N2,H,ALPHA,BETA
C
	NN2 = N
	N21 = N1
	N22 = N2
	
       WRITE(10) NN2,N21,N22
       

       ALPH = INT( ALPHA * 60.)
       BET  = INT( BETA  * 60.)
       H2   = INT( H     * 32768.)
       T2   = INT( X0(1,1) * 32768. / 10.)
*
*	NOTE:  THE LIMITING VALUES  WILL BE - ALPHA, BETA => 1 MINUTE OF ARC
*                                                       H <= 1.
*	                                               T2 <= 10.
*			 DISTANCES AND VELOCITIES <=  10
*					       CP AND CM  <= 10.
	WRITE(10) ALPH, BET, H2, T2
*       WRITE(10) ALPHA,BETA,H,X0(1,1)
    
C
       DO 10 I=1,N
       	DO J = 1, 6
	 XZZ = ( X0(I,J+1) * 32768. / 10.)
	 IF ((XZZ .LT. 32767) .AND. (XZZ.GT.-32767) ) THEN
	  X(J) = INT( XZZ )
	 ELSE
	  IF (XZZ .GT. 32767) THEN
	   X(J) =  32767
	  ELSE
	   X(J) = -32767
	  ENDIF
	 ENDIF
        ENDDO

        WRITE(10) X(1),X(3),X(5),X(2),X(4),X(6)
10     CONTINUE

	DO J=1,6
	 CP2(J) = INT( CP(J) * 32768. / 10.)
	 CM2(J) = INT( CM(J) * 32768. / 10.)
	ENDDO
       WRITE(10) CP2(1),CP2(2),CP2(3),CP2(4),CP2(5),CP2(6)
       WRITE(10) CM2(1),CM2(2),CM2(3),CM2(4),CM2(5),CM2(6)
C
C

        RETURN
        END



C-----------------------------------------------------------------
C
       SUBROUTINE RK4(H)
       DIMENSION XE(10001,7),X0(10001,7),CP(6),CM(6),
     + X(7),F(7),XN(7),XN1(4,7)
       REAL*8 X0,XE,X,CP,CM,H,F,XN,XN1
       INTEGER I,J,N
       COMMON /POS/ X0,XE,CM,CP,N
C
C
C       THIS IS A TRICKY THING- FIRST SOLVE
C       THE MOTION OF THE MASSIVE PARTICLE FOR
C       THE FOUR CALLS, THEN SOLVE FOR EVERYTHING
C       ELSE.
C
C       THE PROBLEM EXISTS BECAUSE WE ARE SOLVING A
C       SET OF SIMULTANIOUS DIFF. EQ. WHICH ARE LINKED
C       TOGETHER.  DOING EACH PARTICLE SEPARATELY WILL
C       NOT WORK, SINCE THE PERTURBING MASS IS CHANGING
C       POSITIONS AT EACH STEP.
C
       DO 100 J=1,7
       X(J) = X0(N,J)
       XN1(1,J)=X0(N,J)
100    XN(J) =X0(N,J)
C
        CALL DIFFEQ(F,X,XN,N)
        DO 101 J=1,7
        XE(N,J)=X0(N,J)+H*F(J)/6.0
        X(J)=X0(N,J)+H*F(J)/2.0
        XN(J)=X(J)
        XN1(2,J)=XN(J)
101     CONTINUE
C
        CALL DIFFEQ(F,X,XN,N)
        DO 102 J=1,7
        XE(N,J)=XE(N,J)+H*F(J)/3.0
        X(J)=X0(N,J)+H*F(J)/2.0
        XN(J)=X(J)
        XN1(3,J)=XN(J)
102     CONTINUE
C
        CALL DIFFEQ(F,X,XN,N)
        DO 103 J=1,7
        XE(N,J)=XE(N,J)+H*F(J)/3.0
        X(J)=X0(N,J)+H*F(J)
        XN(J)=X(J)
        XN1(4,J)=XN(J)
103     CONTINUE
C
        CALL DIFFEQ(F,X,XN,N)
        DO 104 J=1,7
        XE(N,J)=XE(N,J)+H*F(J)/6.0
104     CONTINUE
C
C
       DO 120 I=1,N
       DO 20 J=1,7
        XN(J)=XN1(1,J)
20      X(J)=X0(I,J)
C
       CALL DIFFEQ(F,X,XN,I)
       DO 30 J=1,7
        XE(I,J)=X0(I,J)+H*F(J)/6.0
        XN(J)=XN1(2,J)
30      X(J)=X0(I,J)+H*F(J)/2.0
C
       CALL DIFFEQ(F,X,XN,I)
        DO 40 J=1,7
         XE(I,J)=XE(I,J)+H*F(J)/3.0
         XN(J)=XN1(3,J)
40       X(J)=X0(I,J)+H*F(J)/2.0
C
       CALL DIFFEQ(F,X,XN,I)
         DO 50 J=1,7
          XE(I,J)=XE(I,J)+H*F(J)/3.0
          XN(J)=XN1(4,J)
50        X(J)=X0(I,J)+H*F(J)
C
       CALL DIFFEQ(F,X,XN,I)
         DO 60 J=1,7
60       XE(I,J)=XE(I,J)+H*F(J)/6.0
C
120     CONTINUE
C
         RETURN
         END
C
C
C----------------------------------------------------------------
       SUBROUTINE DIFFEQ(F,X,XN,I)
       DIMENSION F(7),X(7),CP(6),CM(6),XN(7)
       REAL*8 CP,CM,F,X,R1,R21,R2,R22,M1,M2,A1,A2,A3,XN
       REAL*8 X0(10001,7),XE(10001,7)
       INTEGER I,N,K
       COMMON /POS/ X0,XE,CM,CP,N
C
C
       R22=(X(2)-XN(2))*(X(2)-XN(2))+(X(4)-XN(4))*
     1 (X(4)-XN(4))+(X(6)-XN(6))*(X(6)-XN(6))
       R21=X(2)*X(2)+X(4)*X(4)+X(6)*X(6)
       R2N=XN(2)*XN(2)+XN(4)*XN(4)+XN(6)*XN(6)
C
C
C
       R2=SQRT(R22)
       R1=SQRT(R21)
       RN=SQRT(R2N)
C
C
       M1=CM(1)
       M2=CP(1)
C
C      FIND THE MASSES FOR THE POTENTIALS.  IN THIS SIMULATION,
C      CM(1),CP(1) ARE THE MASSES INSIDE A CAUSTIC, CP(2),CM(2)
C      ARE MASSES OF THE CAUSTIC.  THE CAUSTICS ARE AT RADIUS
C      CP(3),CM(3).  CP(4),CM(4) ARE THE SOFTENING LENGTH OF THE
C      POTENTIALS.
C
C
c	print*,m1,m2,r21,r22
c	print*,cp(4), cm(4)
       A1=-M1/(R21+CM(4)*CM(4))
C
C
40     A2=-M2/(R22+CP(4)*CP(4))
C
C
C        REMOVE ACCELERATION FOR PERTURBING GALAXY-
C        FROM PERTURBING GALAXY
C
        IF (I.EQ.N) THEN
        A2=0.0
        R2=1.0
        ENDIF
C
        A3=-M2/(R2N+CP(4)*CP(4))
C
C
       F(1)=1.0
       F(2)=X(3)
       F(3)=A1*X(2)/R1+A2*(X(2)-XN(2))/R2+A3*XN(2)/RN
       F(4)=X(5)
       F(5)=A1*X(4)/R1+A2*(X(4)-XN(4))/R2+A3*XN(4)/RN
       F(6)=X(7)
       F(7)=A1*X(6)/R1+A2*(X(6)-XN(6))/R2+A3*XN(6)/RN
C
C
50       RETURN
         END
*
**************************************************************