C RING C C C c astr306 - hw6 - trying to form a ring galaxy C Steven Janowiecki - 2-Dec-2007 C program ring C some variables... character(len=20) :: ingal real wgal(6),wsat(6),r0,v0,cutr,satgal,r integer outrate, i, indexsave, iname2 real w(6), wdot(6), wnew(6), thistime character(len=20) :: fname character(len=20) :: filename character(len=10) :: iname character*10 nstring parameter(filename='filexxx') real wgaldot(6),wgalnew(6),gw(6) real wsatdot(6),wsatnew(6),sw(6) real rgal,rgal2,gd,rsat,rsat2,sd,GM_gal,GM_sat integer pp C INPUT PARAMETERS C //------------------------------------------------\\ C || || C || How many particles? || parameter (npart=10000) C || sets them up || real x(npart), y(npart), z(npart) real vx(npart), vy(npart), vz(npart) C || || C || Starting time is 0. Stopping time? parameter (timestop=10) C || Step size || parameter (timestep=0.001) C || Total number of steps (timestop/timestep) || parameter (totalit=10000) C || || real tsave(totalit), wsave(6,totalit) C || || C || write an output file every how many steps? || outrate = 10 C || || C || set r0 and v0 (for potential/exponential) || r0 = 0.3 v0 = 2.0**0.5 C || || C || cutoff radius for log/kep potential? || set cutr = 3.0 C || || C || where is initial galaxy file? (x,y,z,vx,vy,vz) || ingal = 'mygalaxy.dat' C || || C || set satellite/galaxy mass ratio || satgal = 0.02 C || || C || set initial w for galaxy (x,y,z,vx,vy,vz) || wgal(1)=0 wgal(2)=0 wgal(3)=0 wgal(4)=0 wgal(5)=0 wgal(6)=0 C || || C || set initial w for satellite (x,y,z,vx,vy,vz) || wsat(1)=0.1 wsat(2)=0.1 wsat(3)=2*cutr wsat(4)=0 wsat(5)=0 wsat(6)=-sqrt(2*v0**2*log((2*cutr/r0)**2)) C || ^^this is set to the escape velocity || C || || C \\------------------------------------------------// C END OF INPUT PARAMETERS C string for file-name generation nstring='0123456789' C read in the intial data file of positions and velocities open(19, FILE=ingal,status='OLD') do i = 1, npart read (19,*) x(i), y(i), z(i), vx(i), vy(i), vz(i) enddo close(19) C work in terms of GM rather than G and M GM_gal = ((v0**2)*cutr) GM_sat = GM_gal*satgal C index variable for output files to be saved indexsave=1 C step over as many steps as we need do thistime = 0, timestop, timestep write(*,*) 'starting on',thistime,'/',timestop C write this timestep to file - first construct filename iname(1:1)=nstring(1+indexsave/1000:1+indexsave/1000) iname2=1+MOD(indexsave,1000)/100 iname(2:2)=nstring(iname2:iname2) iname2=1+MOD(indexsave,100)/10 iname(3:3)=nstring(iname2:iname2) iname3=1+MOD(indexsave,10) iname(4:4)=nstring(iname3:iname3) fname=filename(1:4)//iname(1:4) C only write a file 1/outrate of the time r = mod(indexsave,outrate) if (r .eq. 0) then open(16, FILE=fname) write (16,90) wgal(1),wgal(2),wgal(3), + wgal(4),wgal(5),wgal(6) 90 format (7(e12.5,2x)) write (16,95) wsat(1),wsat(2),wsat(3), + wsat(4),wsat(5),wsat(6) 95 format (7(e12.5,2x)) do i=1,npart write (16,100) x(i), y(i), z(i), vx(i), vy(i), vz(i) enddo 100 format (7(e12.5,2x)) close(16) endif C advance galaxy center and satellite center C main galaxy only feels satellite. call derivs(0.0,wgal,wgaldot,2,wgal,wsat,GM_gal,GM_sat,cutr,v0) call RK4(wgal,wgaldot,6,thistime,timestep,wgalnew,2, + wgal,wsat,GM_gal,GM_sat,cutr,v0) C satellite galaxy only feels main galaxy. call derivs(0.0,wsat,wsatdot,1,wgal,wsat,GM_gal,GM_sat,cutr,v0) call RK4(wsat,wsatdot,6,thistime,timestep,wsatnew,1, + wgal,wsat,GM_gal,GM_sat,cutr,v0) C replace old values wgal(1)=wgalnew(1) wgal(2)=wgalnew(2) wgal(3)=wgalnew(3) wgal(4)=wgalnew(4) wgal(5)=wgalnew(5) wgal(6)=wgalnew(6) wsat(1)=wsatnew(1) wsat(2)=wsatnew(2) wsat(3)=wsatnew(3) wsat(4)=wsatnew(4) wsat(5)=wsatnew(5) wsat(6)=wsatnew(6) C loop over every particle do i = 1,npart C construct w for this particle w(1)=x(i) w(2)=y(i) w(3)=z(i) w(4)=vx(i) w(5)=vy(i) w(6)=vz(i) c calculate dydx (wdot) to send to rk4 call DERIVS(0.0,w,wdot,3,wgal,wsat,GM_gal,GM_sat,cutr,v0) C call RK4 to integrate - particles feel sat and gal call RK4(w,wdot,6,thistime,timestep,wnew,3, + wgal,wsat,GM_gal,GM_sat,cutr,v0) c replace old values w(1)=wnew(1) w(2)=wnew(2) w(3)=wnew(3) w(4)=wnew(4) w(5)=wnew(5) w(6)=wnew(6) C store in overall arrays x(i)=w(1) y(i)=w(2) z(i)=w(3) vx(i)=w(4) vy(i)=w(5) vz(i)=w(6) enddo C advance file-saving index number.. indexsave=indexsave+1 enddo C it's over! all steps are finished stop end SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,pp,gw,sw,GM_gal,GM_sat,cutr,v0) PARAMETER (NMAX=10) DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX),DYT(NMAX),DYM(NMAX) integer pp HH=H*0.5 H6=H/6. XH=X+HH DO 11 I=1,N YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE CALL DERIVS(XH,YT,DYT,pp,gw,sw,GM_gal,GM_sat,cutr,v0) DO 12 I=1,N YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE CALL DERIVS(XH,YT,DYM,pp,gw,sw,GM_gal,GM_sat,cutr,v0) DO 13 I=1,N YT(I)=Y(I)+H*DYM(I) DYM(I)=DYT(I)+DYM(I) 13 CONTINUE CALL DERIVS(X+H,YT,DYT,pp,gw,ww,GM_gal,GM_sat,cutr,v0) DO 14 I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE RETURN END subroutine derivs(x,y,dydx,pp,gw,sw,GM_gal,GM_sat,cutr,v0) c returns derivatives dydx at x c w is coming in as (x y z vx vy vz) c and needs to go out as (vx vy vz ax ay az) C C pp - which potentials? 1 - gal only, 2 - sat only, 3 -gal and sat C gw - galaxy omega, current position to use C sw - satellite omega, current position to use C real x,GM_gal,GM_sat,cutr integer pp dimension y(6), dydx(6), gw(6), sw(6) real galx, galy, galz, galr, satx, saty, satz, satr, v0 C the easy part: d/dt(x) = vx, etc dydx(1)=y(4) dydx(2)=y(5) dydx(3)=y(6) C calculating distances... galx = y(1)-gw(1) galy = y(2)-gw(2) galz = y(3)-gw(3) galr = ((galx**2)+(galy**2)+(galz**2))**0.5 satx = y(1)-sw(1) saty = y(2)-sw(2) satz = y(3)-sw(3) satr = ((satx**2)+(saty**2)+(satz**2))**0.5 C figure out which potential(s) to use, and do so!: if (pp .eq. 1) then C main galaxy only! C check which potential to use if (galr .lt. cutr) then C log potential - main galaxy dydx(4)=-2*(v0**2)*galx/(galr**2) dydx(5)=-2*(v0**2)*galy/(galr**2) dydx(6)=-2*(v0**2)*galz/(galr**2) else C keplerian potential dydx(4)=-GM_gal*galx/(galr**3) dydx(5)=-GM_gal*galy/(galr**3) dydx(6)=-GM_gal*galz/(galr**3) endif elseif (pp .eq. 2) then C satellite only - keplerian dydx(4)=-GM_sat*satx/(satr**3) dydx(5)=-GM_sat*saty/(satr**3) dydx(6)=-GM_sat*satz/(satr**3) else C use both galaxy and satellite, check which galaxy potential if (galr .lt. cutr) then dydx(4)=(-2*(v0**2)*galx/(galr**2)) -(GM_sat*satx/(satr**3)) dydx(5)=(-2*(v0**2)*galy/(galr**2)) -(GM_sat*saty/(satr**3)) dydx(6)=(-2*(v0**2)*galz/(galr**2)) -(GM_sat*satz/(satr**3)) else C keplerian potential for gal, also kep from sat dydx(4)=(-GM_gal*galx/(galr**3)) -(GM_sat*satx/(satr**3)) dydx(5)=(-GM_gal*galy/(galr**3)) -(GM_sat*saty/(satr**3)) dydx(6)=(-GM_gal*galz/(galr**3)) -(GM_sat*satz/(satr**3)) endif endif return end