;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; PROCEDURE Galaxy                   ;;;
;;; PROGRAMMERS EXTRAORDINAIRE         ;;;
;;; Darrin Casebeer &                  ;;;
;;; Kevin Tubbs                        ;;;
;;; Date: 6/01/98-8/08/98              ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;This program is a galaxy collision program adapted from the 
;basic code listed in Astronomy magazine, Dec 1988, by M. C. 
;Schroeder and Neil F. Comins. This program allow two galaxies
;to interact via their gravitational fields. Then nuclei of the 
;TARGET galaxy and the INTRUDER galaxy are point masses, and they 
;move according to their gravitational attraction. The TARGET galaxy 
;has a disk of stars. The stars respond to the gravity of the nuclei
;but do not influence their motions.
;
;The stars are distributed in the target galaxy(only) in rings. 
;The total number of stars is the product of the # of rings and
;the # of stars per ring.

PRO GALAXY,M1,M2,NRR,NRS,VX2,VY2,VZ2,NTSPR
LOADCT,13
X=DBLARR(1000)
Y=DBLARR(1000)
VX=DBLARR(1000)
VY=DBLARR(1000)
Z=DBLARR(1000)
VZ=DBLARR(1000)
X2=0
Y2=0
Z2=35
NS=NRR*NRS
DR=20/(NRR-1)
I=0
;Intialize TARGET galaxy mass, position, and velocity.
X1=0
Y1=0
VX1=0
VY1=0
Z1=0
VZ1=0
SF2=2

;Scale the INTRUDER galaxy mass, position, and velocity.
TIME1=0
M2=M2*M1
X2=X2+X1
Y2=Y2+Y1
Z2=Z2+Z1

;Set up intial load, place stars in concentric rings from r=10 to r=30 
;at intervals of dr.
FOR IR = 1,NRR DO BEGIN
  R=10+(IR-1)*DR
  V=SQRT(M1/R)
  TH=(.5*V/R)*(180/!Pi)
  IF R EQ 10 THEN V=0.9*V
  FOR IT=1,NRS DO BEGIN
    T=(IT-1)*360/NRS
    T1=!Pi*(T-TH)/180
    I=I+1
    
    ;Intialize star positions.
    X(I)=R*COS(T/57.2958)+150
    Y(I)=R*SIN(T/57.2958)+100
    VZ(I)=0
    Z(I)=0
    
    ;Intialize star velocities to place them in circular orbits
    ;around the target galaxy.
    VX(I)=-V*SIN(T1)
    VY(I)=V*COS(T1)
  ENDFOR
ENDFOR

CALL_PROCEDURE,"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=1,1 DO BEGIN
    FOR I=1,NS  DO BEGIN
      
      ;Determine the force on a star due to the galactic centers. SF2 below is
      ;used to prevent overflows during force calculations.
      R1=M1/((X(I)-X1)^2+(Y(I)-Y1)^2+(Z(I)-Z1)^2+SF2)^1.5
      R2=M2/((X(I)-X2)^2+(Y(I)-Y2)^2+(Z(I)-Z2)^2+SF2)^1.5
      
     ;Calculate the stars x,y,z accelerations.
      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))
      
      ;Update stars positions and velocities using a time-centered
      ;leap-frog algorithm.
      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
    
    ;Update positions and velocities of the galactic centers.
    S=((X1-X2)^2+(Y1-Y2)^2+(Z1-Z2)^2+SF2)^1.5
    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+1

  ENDFOR
  CALL_PROCEDURE,"MAINBODY",X,Y,Z,X1,Y1,Z1,SF2,X2,Y2,Z2,VX,VY,VZ,AX,AY,AZ,M1,TIME1,M2,NS,I
  XINTERANIMATE, FRAME=TIME1-1, WINDOW=!D.WINDOW & $
ENDFOR
XINTERANIMATE
END