;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; PROCEDURE  GALAXY                  ;;;
;;; PROGRAMMERS EXTRAORDINAIRE         ;;;
;;; Darrin Casebeer &                  ;;;
;;; Kevin Tubbs                        ;;;
;;; Date: 6/01/98-8/08/98              ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

PRO GALAXY,M1,M2,NRR,NRS,X2,Y2,Z2,VX2,VY2,VZ2,NTSPR
I=0L
X=DBLARR(100000L)
Y=DBLARR(100000L)
VX=DBLARR(100000L)
VY=DBLARR(100000L)
Z=DBLARR(100000L)
VZ=DBLARR(100000L)
NS=NRR*NRS
DR=20d/(NRR-1)
X1=0d
Y1=0d
VX1=0d
VY1=0d
Z1=0d
VZ1=0d
SF2=2d
TIME1=0d
M2=M2*M1
X2=X2+X1
Y2=Y2+Y1
Z2=Z2+Z1

; increment rings through number of rings

FOR IR = 1,NRR DO BEGIN
  R=(10d)+(IR-1d)*DR
  V=SQRT(M1/R)
  TH=(.5d*V/R)*(180d/!Pi)
  IF R EQ 10d THEN V=0.9d*V

  ; set up initial position of stars

  FOR IT=1,NRS DO BEGIN
    T=(IT-1)*360d/NRS
    T1=!Pi*(T-TH)/180d
    I=I+1L
    X(I)=R*COS(T/57.2958d)
    Y(I)=R*SIN(T/57.2958d)
    VZ(I)=0d
    Z(I)=0d
    VX(I)=-V*SIN(T1)
    VY(I)=V*COS(T1)

 ENDFOR
ENDFOR

MAINBODY,X,Y,Z,X1,Y1,Z1,SF2,X2,Y2,Z2,VX,VY,VZ,AX,AY,AZ,M1,TIME1,M2,NS,I



LOOP:
;XINTERANIMATE, SET=[300,300,NTSPR], /showload
FOR K=1,NTSPR DO BEGIN
  FOR J=1L,1 DO BEGIN
    FOR I=1L,NS  DO BEGIN
      R1=M1/((X(I)-X1)^2+(Y(I)-Y1)^2+(Z(I)-Z1)^2+SF2)^1.5d
      R2=M2/((X(I)-X2)^2+(Y(I)-Y2)^2+(Z(I)-Z2)^2+SF2)^1.5d

      ; calculate acceleration on each star based on center of mass
      ; of either galaxy

      AX=R1*(X1-X(I))+R2*(X2-X(I))
      AY=R1*(Y1-Y(I))+R2*(Y2-Y(I))
      AZ=R1*(Z1-Z(I))+R2*(Z2-Z(I))

      ;calculate new velocities for each star based on v = at

      VX(I)=VX(I)+AX
      VY(I)=VY(I)+AY
      VZ(I)=VZ(I)+AZ
      X(I)=X(I)+VX(I)
      Y(I)=Y(I)+VY(I)
      Z(I)=Z(I)+VZ(I)
    ENDFOR
    S=((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2+SF2)^1.5d
    AX=(X2-X1)/S
    AY=(Y2-Y1)/S
    AZ=(Z2-Z1)/S
    VX1=VX1+M2*AX
    VY1=VY1+M2*AY
    VZ1=VZ1+M2*AZ
    VX2=VX2-M1*AX
    VY2=VY2-M1*AY
    VZ2=VZ2-M1*AZ
    X1=X1+VX1
    Y1=Y1+VY1
    Z1=Z1+VZ1
    X2=X2+VX2
    Y2=Y2+VY2
    Z2=Z2+VZ2
    TIME1=TIME1+1d

  ENDFOR
MAINBODY,X,Y,Z,X1,Y1,Z1,SF2,X2,Y2,Z2,VX,VY,VZ,AX,AY,AZ,M1,TIME1,M2,NS,I
WAIT,0.05
;XINTERANIMATE, FRAME=TIME1-1,WINDOW=!D.WINDOW
ENDFOR
;XINTERANIMATE
LASTPART:
END