home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d9xx / d969 / ace.lha / ACE / ACE-2.0.lha / PRGS.lha / Astronomy / GalColSim.b next >
Text File  |  1994-01-10  |  14KB  |  562 lines
  1. { Galaxy Collision Simulator 
  2.  
  3.  Galaxy Collision Program
  4.  by M.C. Schroeder and Neil F. Comins.
  5.  Appeared in Astronomy magazine, December 1988
  6.  as an IBM GW-BASIC program. 
  7.  
  8.  Adapted to Macintosh Turbo Pascal 1.1 by David J. Benn
  9.  Date: 4th April - 17th April, 2nd & 3rd May 1991
  10.  
  11.  Adapted for ACE Amiga BASIC by David J. Benn
  12.  Date: 27th,28th November 1992,
  13.        7th January 1993,
  14.        11th October 1993,
  15.        18th December 1993
  16.     
  17.  This program allows 2 galaxies to interact via
  18.  their gravitational fields. The TARGET galaxy
  19.  is a disk and the INTRUDER galaxy is treated as
  20.  a point mass. }
  21.  
  22. const     SF2=2
  23. const    arr_size=1000
  24. const    aspect=0.625
  25. const    true=-1&
  26. const    false=0&
  27. const    SIM=1
  28. const    PARAMS=2
  29.  
  30. dim     X(arr_size),Y(arr_size),Z(arr_size)
  31. dim    VX(arr_size),VY(arr_size),VZ(arr_size)
  32.  
  33. longint    NRR,NRS,NS,DR
  34.     
  35. single    M1,M2
  36.   
  37. single    X1,Y1,Z1
  38. single    X2,Y2,Z2
  39. single  VX1,VY1,VZ1
  40. single  VX2,VY2,VZ2
  41.     
  42. longint    NTSPR,sim_time
  43.  
  44. SUB title
  45.   title$ = "*** Galaxy Collision Simulator ***"
  46.   color 2
  47.   locate 2,40-len(title$)\2
  48.   PRINT title$   
  49. END SUB
  50.  
  51. SUB await_CR
  52.   color 3
  53.   PRINT
  54.   PRINT "Press <return> to begin simulation..."
  55.   while inkey$<>chr$(13):wend
  56. END SUB
  57.  
  58. SUB instructions
  59. string ans
  60.   PRINT
  61.   color 3
  62.   locate csrlin,25
  63.   PRINT "Do you want instructions? (y/n)"
  64.   repeat
  65.     ans = ucase$(inkey$)
  66.   until ans="Y" or ans="N"
  67.   if ans="Y" then 
  68.     cls
  69.     PRINT
  70.     title
  71.     PRINT
  72.     color 1
  73.     PRINT "This is an adaptation of a BASIC program which appeared in" 
  74.     PRINT "Astronomy magazine, December 1988."
  75.     PRINT
  76.     PRINT "The purpose of the program is to show the large scale effects of"
  77.     PRINT "one galaxy (here refered to as the INTRUDER galaxy) passing nearby"
  78.     PRINT "or through another (TARGET) galaxy."
  79.     PRINT
  80.     PRINT "This is a task which until recent times has been the province of"
  81.     PRINT "of the supercomputer due to the sheer processing power required to"
  82.     PRINT "accurately model the gravitational effects of two galaxies - each"
  83.     PRINT "containing tens or hundreds of billions of stars - upon one ";
  84.     PRINT "another."
  85.     PRINT
  86.     PRINT "However, if the number of stars in the target galaxy is limited and"
  87.     PRINT "the intruder galaxy is considered to be a point gravitational"
  88.     PRINT "source, the problem becomes manageable on a microcomputer."
  89.     PRINT
  90.     PRINT "In this simulation the target galaxy consists of concentric rings" 
  91.     PRINT "of stars."
  92.     await_CR
  93.     cls
  94.     PRINT
  95.     title
  96.     PRINT
  97.     color 1
  98.     PRINT "Before a simulation begins, you will be asked to enter a number of"
  99.     PRINT "parameters. Here is some background information:"
  100.     PRINT
  101.     PRINT "  - 1 parsec = 3.26 light years = approx 30 million million kms"
  102.     PRINT "  - 1 solar mass = the mass of our sun = 1990 thousand million"
  103.     PRINT "    million million million kilograms"
  104.     PRINT "  - Each time step represents 1.2 million years"
  105.     PRINT "  - Each integer X, Y or Z step represents 500 parsecs"
  106.     PRINT "  - Vx, Vy and Vz represent the velocities in the 3 axes"
  107.     PRINT "  - The mass of the target galaxy equals 20 billion solar masses"
  108.     PRINT "  - The mass of the intruder galaxy is entered as a fraction of"
  109.     PRINT "    the target galaxy (eg: 0.25 x 20 billion = 5 billion)."   
  110.     await_CR
  111.     cls
  112.     PRINT
  113.     title
  114.     PRINT
  115.     color 1
  116.     PRINT "The following are some sample initial conditions and their results:"
  117.     PRINT
  118.     PRINT "Result        Intruder Mass   X   Y   Z   Vx   Vy   Vz   Time Steps"
  119.     PRINT
  120.     PRINT "Ring Galaxy       1           7.5  0  35   0    0   -1     65"
  121.     PRINT "Bridge Galaxy     0.25       40   10  10  -1    0    0    120"
  122.     PRINT "Whirlpool Galaxy  0.25      -30   30   0   0   -0.5  0.5  175"
  123.     PRINT
  124.     PRINT "For best results, use about 10 concentric rings and 20 stars per" 
  125.     PRINT "ring. This shows enough detail without slowing the simulation down"
  126.     PRINT "too much."    
  127.     PRINT 
  128.     PRINT "In what follows, you will also be given the opportunity to choose"
  129.     PRINT "one of the above 3 resultant galaxies rather than entering initial"
  130.     PRINT "condition values." 
  131.     PRINT
  132.     PRINT "This program may also be run from the shell/CLI, with the following"
  133.     PRINT "syntax:"
  134.     PRINT
  135.     PRINT "    GalColSim [ring | bridge | whirlpool]"
  136.     await_CR
  137.   end if
  138. END SUB
  139.  
  140. SUB entry_method
  141. string ans 
  142.   cls
  143.   color 1
  144.   PRINT:PRINT
  145.   PRINT "Enter [P]arameters or [N]ame of simulation?"
  146.   repeat
  147.     ans=ucase$(inkey$)
  148.   until ans="P" or ans="N"
  149.   case
  150.     ans="P" : entry_method=PARAMS
  151.     ans="N" : entry_method=SIM
  152.   end case
  153. END SUB
  154.  
  155. SUB set_parameters(opt)
  156. shared      NRR,NRS,NS,DR
  157. shared     M2
  158. shared     X2,Y2,Z2
  159. shared     VX2,VY2,VZ2
  160. shared     NTSPR
  161.   '..initialise parameters
  162.   if opt=1 then
  163.     '..ring
  164.     M2=1
  165.     X2=7.5 : Y2=0  : Z2=35
  166.     VX2=0  : VY2=0 : VZ2=-1    
  167.     NTSPR=65
  168.   else
  169.     if opt=2 then
  170.       '..bridge
  171.       M2=0.25
  172.       X2=40  : Y2=10 : Z2=10
  173.       VX2=-1 : VY2=0 : VZ2=0
  174.       NTSPR=120   
  175.     else
  176.       '..whirlpool
  177.       M2=0.25
  178.       X2=-30 : Y2=30    : Z2=0
  179.       VX2=0  : VY2=-0.5 : VZ2=0.5
  180.       NTSPR=175
  181.     end if
  182.   end if 
  183.  
  184.   '..number of rings and stars per ring
  185.   NRR=10
  186.   NRS=20
  187.   NS=NRR*NRS
  188.   if NRR<2 then NRR=2
  189.   DR=20 \ (NRR-1)
  190. END SUB
  191.  
  192. SUB args_present
  193. string end_state
  194.   if argcount<>1 then
  195.     args_present=false
  196.     exit sub
  197.   else 
  198.     args_present=true
  199.     end_state=ucase$(arg$(1))
  200.     case
  201.       end_state="RING"      : opt=1
  202.       end_state="BRIDGE"    : opt=2
  203.       end_state="WHIRLPOOL" : opt=3
  204.     end case
  205.     set_parameters(opt)       
  206.   end if
  207. END SUB
  208.  
  209. SUB get_simulation
  210. shortint opt
  211.  
  212.   '..ask for simulation end-state
  213.   PRINT
  214.   PRINT "Which end-state?"
  215.   PRINT
  216.   PRINT "1. Ring Galaxy"
  217.   PRINT "2. Bridge Galaxy"
  218.   PRINT "3. Whirlpool Galaxy"
  219.   repeat 
  220.     opt=val(inkey$)
  221.   until opt>=1 and opt<=3
  222.   
  223.   set_parameters(opt)
  224. END SUB
  225.  
  226. SUB get_parameters
  227. {The stars are distributed in the
  228.  TARGET galaxy in rings. The total number
  229.  of stars is the product of the # of rings and
  230.  the # of stars per ring.}
  231. shared     NRR,NRS,NS,DR
  232. shared    M2
  233. shared    X2,Y2,Z2
  234. shared    VX2,VY2,VZ2
  235. shared    NTSPR
  236. string     ANS  
  237.  repeat
  238.   cls
  239.   PRINT
  240.   PRINT "Enter the initial conditions..."
  241.   PRINT
  242.   PRINT "Input the number of rings of stars in the TARGET galaxy (try 8-20):"
  243.   input NRR
  244.   PRINT
  245.   PRINT "Input the number of stars per ring in the TARGET galaxy (try 15-30):"
  246.   input NRS
  247.   PRINT
  248.   NS=NRR*NRS
  249.   if NRR<2 then NRR=2
  250.   DR=20 \ (NRR-1)
  251.   PRINT "Input the mass fraction of the INTRUDER galaxy"
  252.   PRINT "in units of the TARGET galaxy mass (try 0.25-1):"
  253.   input M2
  254.   PRINT
  255.   PRINT "Input the initial X Y Z co-ordinates of the INTRUDER galaxy"
  256.   PRINT "(eg: 40 10 10):"
  257.   PRINT "The TARGET galaxy is located at 0 0 0."
  258.   PRINT "X: "; 
  259.   input X2
  260.   PRINT "Y: "; 
  261.   input Y2
  262.   PRINT "Z: "; 
  263.   input Z2
  264.   PRINT
  265.   PRINT "Input the initial X Y Z  velocities"
  266.   PRINT "of the INTRUDER galaxy (eg: -1 0 0):"
  267.   PRINT "The TARGET galaxy is initially at rest."
  268.   PRINT "X: "; 
  269.   input VX2
  270.   PRINT "Y: "; 
  271.   input VY2
  272.   PRINT "Z: "; 
  273.   input VZ2
  274.   {not too small (0 is ok) - system may crash at around +- .34 or so...}
  275.   if VX2<0.5 then VX2=VX2*2;
  276.   if VY2<0.5 then VY2=VY2*2;
  277.   if VZ2<0.5 then VZ2=VZ2*2;
  278.   PRINT
  279.   {the position of the TARGET galaxy
  280.    is assigned by the computer}
  281.   PRINT "Input number of time steps for this run (try 50-200)"
  282.   input NTSPR
  283.   PRINT
  284.   PRINT "Are these inputs correct? (y/n) "
  285.   repeat
  286.     ANS = ucase$(inkey$)
  287.   until ANS="Y" or ANS="N"
  288.  until ANS="Y"
  289. END SUB
  290.     
  291. SUB update_display
  292. {update screen display}
  293. shared     X,Y,VX,VY,Z,VZ
  294. shared    NRR,NRS,NS,DR    
  295. shared    M1,M2  
  296. shared    X1,Y1,Z1
  297. shared    X2,Y2,Z2
  298. shared  VX1,VY1,VZ1
  299. shared  VX2,VY2,VZ2   
  300. shared    NTSPR,sim_time
  301. single     XC,YC,ZC,XX1,YY1,ZZ1,XX2,YY2,ZZ2
  302. longint    I,J
  303. single    XP,YP,ZP
  304. string    ch
  305.     {calculate centre of mass of system
  306.      for use as centre of output display}
  307.     XC=(M1*X1+M2*X2)/(M1+M2)
  308.     YC=(M1*Y1+M2*Y2)/(M1+M2)
  309.     ZC=(M1*Z1+M2*Z2)/(M1+M2)
  310.     {set up output display - clear text & graphics}
  311.     cls
  312.     {print out display headings}  
  313.     color 3
  314.     locate 2  
  315.     PRINT "        Polar             Edge-on"
  316.     color 3
  317.     locate 22
  318.     PRINT "        Step:";sim_time
  319.     color 5
  320.     locate 23
  321.     PRINT "   Real Time:";int((sim_time+1)*1.2);"million years"
  322.     color 4
  323.     locate 24
  324.     PRINT "   Press any key to stop"
  325.     '..view frames
  326.     LINE (0,0)-(159,166),2,b
  327.     LINE (161,0)-(319,166),2,b
  328.     {calculate position of galactic centres and display on screen}
  329.     XX1=(X1-XC)
  330.     YY1=(Y1-YC)
  331.     ZZ1=(Z1-ZC)
  332.     XX2=(X2-XC)
  333.     YY2=(Y2-YC)
  334.     ZZ2=(Z2-ZC)
  335.     ZP=Z(1)-ZC + Z(NS)-ZC
  336.     CIRCLE(int(XX1+75),int(YY1+75)),2,1,,,aspect    {target x-y}
  337.     CIRCLE(int(XX1+225),int(ZP+75)),2,1,,,aspect    {target x-z}
  338.     CIRCLE(int(XX2+75),int(YY2+75)),1,5,,,aspect         {intruder x-y}
  339.     CIRCLE(int(XX2+225),int(2*ZZ2+75)),1,5,,,aspect       {intruder x-z}
  340.     {place stars on screen}
  341.     FOR I=1 TO NS
  342.       XP=int(X(I)-XC)
  343.       YP=int(Y(I)-YC)
  344.       ZP=2*(Z(I)-ZC)
  345.       PSET(int(XP+75),int(YP+75)),1
  346.       PSET(int(XP+225),int(ZP+75)),1
  347.     NEXT
  348. END SUB
  349.  
  350. SUB define_galaxies
  351. shared     X,Y,VX,VY,Z,VZ
  352. shared    NRR,NRS,NS,DR    
  353. shared    M1,M2  
  354. shared    X1,Y1,Z1
  355. shared    X2,Y2,Z2
  356. shared  VX1,VY1,VZ1
  357. shared  VX2,VY2,VZ2   
  358. shared    NTSPR,sim_time
  359. longint    IR,IT,I
  360. single    R,RV,T,T1,TH,V
  361.  PRINT
  362.  PRINT "*** Defining structure & position" 
  363.  PRINT "    of galaxies. Please wait."
  364.  {initialise TARGET galaxy mass,
  365.   position and velocity}
  366.  M1=5
  367.  X1=150
  368.  Y1=100
  369.  Z1=0
  370.  VX1=0
  371.  VY1=0
  372.  VZ1=0
  373.  {scale the INTRUDER galaxy mass,
  374.   position and velocities.}
  375.  sim_time=0
  376.  M2=M2*M1
  377.  X2=X2+X1
  378.  Y2=Y2+Y1
  379.  Z2=Z2+Z1
  380.  {setup initial load
  381.   place stars in concentric rings
  382.   from r=10 to r=30 at intervals of DR}
  383.  FOR IR=1 TO NRR
  384.     R=10+(IR-1)*DR
  385.     V=SQR(M1/R)
  386.     TH=(0.5*V/R)*(180/3.14159)
  387.     IF R=10 THEN V=0.9*V
  388.     FOR IT=1 TO NRS
  389.         T=(IT-1)*360/NRS
  390.         T1=3.14159*(T-TH)/180
  391.         I=I+1      
  392.         {initialise star positions}
  393.         X(I)=R*COS(T/57.2958)+150
  394.         Y(I)=R*SIN(T/57.2958)+100
  395.         VZ(I)=0
  396.         Z(I)=0
  397.         {initialise star velocities to place them in circular
  398.          orbits about the TARGET galaxy.}
  399.         VX(I)=-V*SIN(T1)
  400.         VY(I)=V*COS(T1)
  401.     NEXT
  402.  NEXT
  403. END SUB
  404.  
  405. SUB particle_pusher
  406. shared     X,Y,VX,VY,Z,VZ
  407. shared    NRR,NRS,NS,DR    
  408. shared    M1,M2  
  409. shared    X1,Y1,Z1
  410. shared    X2,Y2,Z2
  411. shared  VX1,VY1,VZ1
  412. shared  VX2,VY2,VZ2   
  413. shared    NTSPR,sim_time
  414. longint    I,J,K
  415. single    R1,R2
  416. single    AX,AY,AZ
  417. single    RR 
  418. single    Ex,Ey,Ez
  419.     FOR K=1 TO NTSPR 
  420.         FOR J=1 TO 1 
  421.             FOR I=1 TO NS
  422.              {determine the force on a star
  423.               due to the galactic centres
  424.               SF2 below, called the softening factor,
  425.               is used to prevent overflows
  426.               during force calculations
  427.               SF2 is assigned a value above.
  428.               You may experiment with its value
  429.               but it should be kept as small as possible
  430.               to better approximate a true 1/r**2 force.}
  431.               Ex=(X(I)-X1)^2
  432.               Ey=(Y(i)-Y1)^2
  433.               Ez=(Z(i)-Z1)^2
  434.               R1=M1/(Ex+Ey+Ez+SF2)^1.5
  435.               Ex=(X(I)-X2)^2
  436.               Ey=(Y(i)-Y2)^2
  437.               Ez=(Z(i)-Z2)^2
  438.               R2=M2/(Ex+Ey+Ez+SF2)^1.5
  439.               {calculate stars x,y,z accelerations}
  440.               AX=R1*(X1-X(I))+R2*(X2-X(I))
  441.               AY=R1*(Y1-Y(I))+R2*(Y2-Y(I))
  442.               AZ=R1*(Z1-Z(I))+R2*(Z2-Z(I))
  443.               {update star positions and velocities
  444.                using a time centred leap-frog algorithm}
  445.               VX(I)=VX(I)+AX
  446.               VY(I)=VY(I)+AY
  447.               VZ(I)=VZ(I)+AZ
  448.               X(I)=X(I)+VX(I)
  449.               Y(I)=Y(I)+VY(I)
  450.               Z(I)=Z(I)+VZ(I)
  451.               if inkey$<>"" then exit sub
  452.             NEXT
  453.           {update positions
  454.            and velocities of the galactic centres}
  455.            Ex=(X1-X2)^2
  456.            Ey=(Y1-Y2)^2
  457.            Ez=(Z1-Z2)^2
  458.            RR=(Ex+Ey+Ez+SF2)^1.5
  459.            AX=(X2-X1)/RR
  460.            AY=(Y2-Y1)/RR
  461.            AZ=(Z2-Z1)/RR
  462.            VX1=VX1+M2*AX
  463.            VY1=VY1+M2*AY
  464.            VZ1=VZ1+M2*AZ
  465.            VX2=VX2-M1*AX
  466.            VY2=VY2-M1*AY
  467.            VZ2=VZ2-M1*AZ
  468.            X1=X1+VX1
  469.            Y1=Y1+VY1
  470.            Z1=Z1+VZ1
  471.            X2=X2+VX2
  472.            Y2=Y2+VY2
  473.            Z2=Z2+VZ2
  474.            ++sim_time
  475.         NEXT
  476.         update_display
  477.      NEXT   
  478. END SUB
  479.  
  480. SUB continue
  481. shared     NTSPR
  482. string     ans
  483.     color 2
  484.     locate 24
  485.     PRINT "   Continue simulation? (y/n)"
  486.     repeat
  487.       ans = ucase$(inkey$)
  488.     until ans="Y" or ans="N"
  489.     if ans="Y" then 
  490.       continue=true
  491.       locate 24
  492.       PRINT "   Enter number of time steps: ";
  493.       input NTSPR
  494.     else 
  495.       continue=false
  496.     end if
  497. END SUB
  498.  
  499. SUB finished 
  500. string ans
  501.   cls
  502.   color 5
  503.   locate 12,8
  504.   PRINT "Another Simulation? (y/n)"
  505.   repeat
  506.     ans=ucase$(inkey$)
  507.   until ans="Y" or ans="N"
  508.   if ans="N" then finished=true else finished=false
  509. END SUB
  510.  
  511. SUB open_hires
  512.   SCREEN 2,640,225,2,2
  513.   palette 0,0,0,0    '..black
  514.   palette 1,1,1,1    '..white
  515.   palette 2,1,1,0.13    '..yellow
  516.   palette 3,1,0,0    '..red
  517.   color 1
  518. END SUB
  519.  
  520. SUB close_hires
  521.   SCREEN CLOSE 2
  522. END SUB
  523.  
  524. { ** main ** }  
  525. SCREEN 1,320,225,3,1
  526.  
  527. palette 0,0,0,0        '..black
  528. palette 1,1,1,1        '..white
  529. palette 2,1,0,0        '..red
  530. palette 3,0,1,0        '..green
  531. palette 4,0,0,1        '..blue
  532. palette 5,1,1,0.13    '..yellow
  533.  
  534. color 1
  535.  
  536.     argument = args_present
  537.  
  538.     repeat  
  539.       if not argument then 
  540.         open_hires
  541.         title
  542.         instructions 
  543.         if entry_method = SIM then
  544.           get_simulation
  545.         else
  546.           get_parameters
  547.         end if
  548.         define_galaxies 
  549.         close_hires
  550.         update_display
  551.       else
  552.         define_galaxies
  553.     update_display
  554.         argument=false
  555.       end if
  556.       repeat   
  557.         particle_pusher
  558.       until not continue
  559.     until finished
  560.  
  561. SCREEN CLOSE 1
  562.