AngMon: Two particle system. 4 initial conditions. Angular momentum xl(), xl(4) magnitude. Set 1: 0.6510 0.5850 -0.2380 -0.7550 -0.8280 -0.8650 -0.7260 0.9310 -0.0960 0.0000 0.3570 -0.2090 0.1070 -0.6600 XMM,xmu: 1.582000 0.383111 r12: 0.681000 -0.238000 -1.112000 1.325500 v12: -0.619000 -0.972000 -0.066000 1.154253 vcm: -0.463721 -0.292982 -0.687159 0.879241 xl,E: -0.408073 0.280925 -0.310034 0.584434 -0.202038 Second Calculation (contol check): xl,E: -0.408073 0.280925 -0.310034 0.584434 -0.202038 e,p: 0.138826 1.471011 Set 2: 1.5100 0.4600 -0.3590 -0.2340 -0.9180 -0.9410 -0.3230 0.1260 -0.0660 -0.0900 -0.8090 0.7890 0.7880 0.6200 XMM,xmu: 1.636000 0.116296 r12: 0.526000 -0.269000 0.575000 0.824416 v12: -1.707000 -1.729000 -0.943000 2.606250 vcm: -0.786532 -0.807837 -0.250373 1.154955 xl,E: 0.145119 -0.056462 -0.159167 0.222669 0.164191 Second Calculation (contol check): xl,E: 0.145119 -0.056462 -0.159167 0.222669 0.164191 e,p: 2.206259 2.240828 Set 3: 1.3280 -0.1250 0.8980 0.1940 -0.4520 0.1720 0.1250 1.9990 -0.4490 -0.0850 -0.4540 -0.9760 -0.9900 -0.9680 XMM,xmu: 3.327000 0.797918 r12: 0.324000 0.983000 0.648000 1.221134 v12: 0.524000 1.162000 1.093000 1.679127 vcm: -0.766841 -0.526178 -0.531720 1.071277 xl,E: 0.256485 -0.011634 -0.110595 0.279555 -1.049087 Second Calculation (contol check): xl,E: 0.256485 -0.011634 -0.110595 0.279555 -1.049087 e,p: 0.985312 0.036895 Set 4: 0.1800 0.2040 -0.9680 -0.7530 -0.8110 -0.6320 0.7840 1.5600 -0.8890 -0.9790 0.8540 -0.3230 -0.7740 -0.5330 XMM,xmu: 1.740000 0.161379 r12: 1.093000 0.011000 -1.607000 1.943507 v12: -0.488000 0.142000 1.317000 1.411665 vcm: -0.373483 -0.759310 -0.396759 0.934590 xl,E: 0.039164 -0.105746 0.025913 0.115705 0.016317 Second Calculation (contol check): xl,E: 0.039164 -0.105746 0.025913 0.115705 0.016317 e,p: 1.017022 0.295431 Fortran Program: program AngMom ! Bernd Berg, Feb 12 2011. C Angular Momentum in CM from fro given initial conditions. C LL: Landau-Lifshitz Mechanics. implicit real*8 (a-h,o-z) logical ltest character c*2 C n initial conditions, dimension 3, 4 magnitude, 2 particles. parameter(n=4,ngnu=1000,iuo=6,iud0=10,iud1=11,iud2=12) parameter(zero=1.d00,one=1.d00) dimension r(4,2),v(4,2),xm(2), r12(4),v12(4),rcm(4),vcm(4),xl(4) ltest=.true. ltest=.false. if(ltest) stop "AngMon: ltest." C write(iuo,'(/,"AngMon: Two particle system.")') write(iuo,'(/,I3," initial conditions.")') n write(iuo,'(/," Angular momentum xl(), xl(4) magnitude.")') pi=acos(-one) C Initial values from file: open(iud0,file="Kepler.txt",form='formatted',status='old') open(iud1,file="Orbits0.d",form='formatted',status='unknown') do i=1,n write(iuo,'(/,"Set",I2,":",/)') i read(iud0,*) i_in,xm(1),(r(id,1),id=1,3),(v(id,1),id=1,3) write(iuo,'(7F9.4)') xm(1),(r(id,1),id=1,3),(v(id,1),id=1,3) read(iud0,*) i_in,xm(2),(r(id,2),id=1,3),(v(id,2),id=1,3) write(iuo,'(7F9.4)') xm(2),(r(id,2),id=1,3),(v(id,2),id=1,3) XMM=xm(1)+xm(2) ! Total mass (mu in LL). xmu=xm(1)*xm(2)/XMM ! Reduced mass (m in LL). write(iuo,'(/,4X,"XMM,xmu:",10X,2F13.6)') XMM,xmu do id=1,3 r12(id)=r(id,1)-r(id,2) v12(id)=v(id,1)-v(id,2) vcm(id)=(xm(1)*v(id,1)+xm(2)*v(id,2))/XMM enddo r12(4)=sqrt(r12(1)**2+r12(2)**2+r12(3)**2) v12(4)=sqrt(v12(1)**2+v12(2)**2+v12(3)**2) vcm(4)=sqrt(vcm(1)**2+vcm(2)**2+vcm(3)**2) write(iuo,'(4X,"r12: ",4F13.6)') r12 write(iuo,'(4X,"v12: ",4F13.6)') v12 write(iuo,'(4X,"vcm: ",4F13.6)') vcm xl(1)=xmu*(r12(2)*v12(3)-r12(3)*v12(2)) xl(2)=xmu*(r12(3)*v12(1)-r12(1)*v12(3)) xl(3)=xmu*(r12(1)*v12(2)-r12(2)*v12(1)) xl(4)=sqrt(xl(1)**2+xl(2)**2+xl(3)**2) E=xmu*v12(4)**2/2-xm(1)*xm(2)/r12(4) write(iuo,'(4X,"xl,E:",5F13.6)') xl,E write(iuo,'(4X,"Second Calculation (contol check):")') do id=1,3 ! CM frame coordinates. r(id,1)=r(id,1)-rcm(id) r(id,2)=r(id,2)-rcm(id) v(id,1)=v(id,1)-vcm(id) v(id,2)=v(id,2)-vcm(id) enddo r(4,1)=sqrt(r(1,1)**2+r(2,1)**2+r(3,1)**2) r(4,2)=sqrt(r(1,2)**2+r(2,2)**2+r(3,2)**2) v(4,1)=sqrt(v(1,1)**2+v(2,1)**2+v(3,1)**2) v(4,2)=sqrt(v(1,2)**2+v(2,2)**2+v(3,2)**2) xl(1)=xm(1)*(r(2,1)*v(3,1)-r(3,1)*v(2,1)) & +xm(2)*(r(2,2)*v(3,2)-r(3,2)*v(2,2)) xl(2)=xm(1)*(r(3,1)*v(1,1)-r(1,1)*v(3,1)) & +xm(2)*(r(3,2)*v(1,2)-r(1,2)*v(3,2)) xl(3)=xm(1)*(r(1,1)*v(2,1)-r(2,1)*v(1,1)) & +xm(2)*(r(1,2)*v(2,2)-r(2,2)*v(1,2)) E=xm(1)*v(4,1)**2/2+xm(2)*v(4,2)**2/2-xm(1)*xm(2)/r12(4) write(iuo,'(4X,"xl,E:",5F13.6)') xl,E C Prepare plots of effective potentials: write(c,'(I2.2)') i open(iud2,file="Ueff"//c//".d", & form='formatted',status='unknown') do ignu=1,ngnu rr=0.01d00*ignu Ueff=-xm(1)*xm(2)/rr+(xl(4)/rr)**2/(2*xmu) write(iud2,'(3F13.6)') rr,Ueff,E enddo close(iud2) C Prepare orbit plots: alpha=xm(1)*xm(2) e=sqrt(one+2*E*xl(4)**2/alpha**2/xmu) ! Eccentricity. p=xl(4)**2/alpha/xmu ! Half latus rectum. theta0=acos((p/r12(4)-one)/e) ! Initial conditions. rdv=r12(1)*v12(1)+r12(2)*v12(2)+r12(3)*v12(3) if(rdv.lt.zero) theta0=-theta0 x0=r12(4)*cos(theta0) y0=r12(4)*sin(theta0) thetamax=pi dtheta=2*thetamax/(2*ngnu) if(e.ge.one) then thetamax=acos(-one/e) dtheta=2*thetamax/(2*ngnu+1) endif open(iud2,file="Orbit"//c//".d", & form='formatted',status='unknown') do ignu=-ngnu,ngnu theta=ignu*dtheta rr=p/(one+e*cos(theta)) x=rr*cos(theta) y=rr*sin(theta) write(iud2,'(4F13.6)') x0,y0,x,y enddo close(iud2) C enddo close(iud0) close(iud1) stop end