c astr306 hw4 orbit nbody (1body) C Steven Janowiecki - 11/9/2007 C C revised for HW6 - from the ground up... C C (good version!! don't use myorbit.f) C C C UNITS ARE A MAJOR PROBLEM STILLLLLLLLLL C C program myorbit2 real x(10000), y(10000), z(10000) real vx(10000), vy(10000), vz(10000) integer i, indexsave parameter (totalit=10000) parameter (timestep=0.001) parameter (timestop=10) real tsave(totalit) real w(6), wdot(6), wnew(6), thistime, wsave(6,totalit) character(len=20) :: fname character(len=20) :: filename character(len=10) :: iname integer iname2 character*10 nstring parameter(filename='filexxx') real r real wgal(6),wgaldot(6),wgalnew(6),gw(6) real wsat(6),wsatdot(6),wsatnew(6),sw(6) real rgal,rgal2,gd,rsat,rsat2,sd real GM_gal,GM_sat real cutr real r0 parameter(cutr = 3.0) integer pp nstring='0123456789' C read in the intial data file of positions and velocities open(19, FILE='mygalaxy.dat',status='OLD') do i = 1, 10000 read (19,*) x(i), y(i), z(i), vx(i), vy(i), vz(i) enddo close(19) C from exponential disk r0 = 0.3 C v0 ^2 * R GM_gal = (((2.0**0.5)**2)*cutr) GM_sat = GM_gal*0.02 C initial w for galaxy and satellite wgal(1)=0 wgal(2)=0 wgal(3)=0 wgal(4)=0 wgal(5)=0 wgal(6)=0 wsat(1)=0.1 wsat(2)=0.1 wsat(3)=6.0 wsat(4)=0 wsat(5)=0 C v_esc = sqrt( 2 * phi(2 * cutr)) C v_esc = sqrt(2 * (sqrt(2))**2 * ln( (10r0/r0)**2)) C v_esc = sqrt(2*2*ln(100))=4.2919 wsat(6)=-4.2919 indexsave=1 C step over as many steps as we need do thistime = 0, timestop, timestep write(*,*) 'starting on',thistime,'/',timestop C write this to file 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 1/10 of the time r = mod(indexsave,10) C r = 0 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,10000 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. start with wgaldot: call derivs(0.0,wgal,wgaldot,2,wgal,wsat,GM_gal,GM_sat,cutr) call RK4(wgal,wgaldot,6,thistime,timestep,wgalnew,2, + wgal,wsat,GM_gal,GM_sat,cutr) C satellite galaxy only feels main galaxy. start with wsatdot: call derivs(0.0,wsat,wsatdot,1,wgal,wsat,GM_gal,GM_sat,cutr) call RK4(wsat,wsatdot,6,thistime,timestep,wsatnew,1, + wgal,wsat,GM_gal,GM_sat,cutr) 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,10000 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) for rk4 call DERIVS(0.0,w,wdot,3,wgal,wsat,GM_gal,GM_sat,cutr) C call RK4 to integrate call RK4(w,wdot,6,thistime,timestep,wnew,3, + wgal,wsat,GM_gal,GM_sat,cutr) 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) 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 indexsave=indexsave+1 enddo stop end SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,pp,gw,sw,GM_gal,GM_sat,cutr) 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) 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) 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) DO 14 I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE RETURN END C X subroutine derivs(x,y,dydx,pp,gw,sw,GM_gal,GM_sat,cutr) 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 real v0 C km/s v0 = (2.0**0.5) dydx(1)=y(4) dydx(2)=y(5) dydx(3)=y(6) 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 write (*,*) y(1),y(2),y(3),satx,saty,satz C write (*,*) y(3),sw(3),satz C figure out which potentials to include: C main galaxy Only 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 C _ _ _ _ _ _ _ _ C _-(_)- _-(_)- _-(_)- _-(_)- _-(_)- _-(_)- _-(_)- _-(_)- C`(___) `(___) `(___) `(___) `(___) `(___) `(___) `(___) C // \\ // \\ // \\ // \\ // \\ // \\ // \\ // \\ C return end