DEFDBL a-z DECLARE FUNCTION r1 (x, y, mu) DECLARE FUNCTION r2 (x, y, mu) DECLARE FUNCTION vxy (x, y, mu) DECLARE FUNCTION dvdx (x, y, mu) DECLARE FUNCTION dvdy (x, y, mu) DECLARE SUB GRAPHSET(mu,xlps(),ylps()) DECLARE SUB Nbranch(mu,x,y,k,dx,dy,Nb%,xnxs(),ynxs()) DECLARE SUB LAGPS(mu,xlps(),xlps()) DIM NCS%(2) DIM mus(2),dxs(2), dys(2),x0s(20,2),y0s(20,2),k0s(20,2),xnxs(5),ynxs(5) DIM xlps(4),ylps(4) ' PRINT CURDIR$ ' SLEEP 'OPEN "C:\FreeBASIC017b\P-Hellings\Chapt5P.dat" FOR INPUT AS #1 'OPEN "C:\FreeBASIC017b\P-Hellings\Chapt5P.out" FOR OUTPUT AS #2 CLS PRINT " P.Hellings: Astrophysics with a PC : EQUIPOTENTIAL CURVES" PRINT " ---------------------------------------------" ' Adapted by Kiyoshi Kawabata from the program of Chapter 5 ' on May 12, 2008 INPUT " 0 for manual data input, 1 for Fig. 5.7, and 2 for Fir. 5.8 : "; Nkey% PRINT "" IF Nkey%=0 THEN PRINT " Input of initial conditions and parameters : " INPUT " mass parameter mu"; mus(0) INPUT " stepsize dx"; dxs(0) INPUT " stepsize dy"; dys(0) NCASES%=0 FOR I%=1 TO 20 INPUT " x0 (Input 100 to quit"; x0s(I%,0) IF(x0s(I%,0)> 99.9999) THEN GOTO 3 INPUT " y0 "; y0s(I%,0) k0s(I%,0)=vxy(x0,yo,mu) NCASES%=NCASES%+1 PRINT "(";NCASES%;")" NEXT I% 3: NCS%(0)=NCASES% ELSE FOR I%=1 TO 2 READ mus(I%),dxs(I%),dys(I%) READ NCS%(I%) FOR J%=1 TO NCS%(I%) READ x0s(J%,I%), y0s(J%,I%) NEXT J% NEXT I% ENDIF NCASES%=NCS%(Nkey%) mu=mus(Nkey%) dxin0=dxs(Nkey%) dyin0=dys(Nkey%) CALL LAGPS(mu,xlps(),ylps()) PRINT " Lagrangian points:" FOR I%=1 TO 4 PRINT " L";I%;" : x=";xlps(I%);" y=";ylps(I%) NEXT I% FOR KCASE%=1 TO NCASES% x0=x0s(KCASE%,Nkey%) y0=y0s(KCASE%,Nkey%) k0 = vxy(x0, y0, mu) k0s(KCASE%,Nkey%)=k0' ' Find the number of branchings at the given initial position for an EP curve: CALL Nbranch(mu,x0,y0,k0,dxin0,dyin0,Nb%,xnxs(),ynxs()) ' Initialize a graphic screen IF KCASE%=1 THEN CALL GRAPHSET(mu,xlps(),ylps()) LOCATE(1,2) PRINT "mu=";mu ENDIF ' Start plottting the KCASE%-th EP curve PSET(x0,y0) FOR Lb%=1 TO Nb% x1=xnxs(Lb%) y1=ynxs(Lb%) LINE (x0,y0)-(x1,y1) dx=x1-x0 dy=y1-y0 ISIGX=SGN(dx) ISIGY=SGN(dy) IF ABS(dx)< 1.e-6 THEN dx=ISIGX%*1.e-6 IF ABS(dy)< 1.e-6 THEN dy=ISIGY%*1.e-6 drc%=1 IF ABS(dy/dx)>1.0 THEN drc%=2 FOR KDUMMY%=1 TO 1000 ' use midpoint method to compute y as function of x if IDRC is 1 IF drc%=1 THEN x12=x1+.5*dx y12=y1-.5*dx*dvdx(x1,y1,mu)/dvdy(x1,y1,mu) x=x1+dx y=y1-dx*dvdx(x12,y12,mu)/dvdy(x12,y12,mu) ' add one corretion of the computed value of y y=y-(vxy(x,y,mu)-k0)/dvdy(x,y,mu) ' use midpoint method to compute x as tion of y if IDRC is 2 ELSE y12=y1+.5*dy x12=x1-.5*dy*dvdy(x1,y1,mu)/dvdx(x1,y1,mu) y=y1+dy x=x1-dy*dvdy(x12,y12,mu)/dvdx(x12,y12,mu) ' add one corretion of the computed value of x x=x-(vxy(x,y,mu)-k0)/dvdx(x,y,mu) ENDIF dxstep=x-x1 dystep=y-y1 IF y<0.0 GOTO 101 ' draw a line from the currrent position to the point (x,y) LINE -(x,y) __SLEEP 1 IF KDUMMY%>10 AND (x-x0)^2+(y-y0)^2 < 1.e-4 GOTO 101 ISIGX%=SGN(dxstep) ISIGY%=SGN(dystep) IF ABS(dxstep)< 1.e-6 THEN dxstep=ISIG%*1.e-6 IF ABS(dystep)< 1.e-6 THEN dystep=ISIG%*1.e-6 drc%=1 IF ABS(dystep/dxstep)>1.0 THEN drc%=2 dx=x-x1 dy=y-y1 x1=x y1=y 100: NEXT KDUMMY% 101: NEXT Lb% LOCATE (2,5*Kcase%) PRINT "(";KCASE%;") NEXT KCASE% COLOR 14 PRINT "press any key to terminate this job" DO LOOP WHILE INKEY$ = "" ' Input data to reproduce Figs. 5.7 and 5.8 ' Fig. 5.7 DATA 0.4, 0.005, 0.005 ': mu, dx, dy DATA 8 ': NEP%, the number of EP curves DATA -1.23081, 0. ': x0, y0 DATA -1.08190, 0. DATA -0.4407780, 0. DATA -0.14162, 0. DATA 0.200, 0. DATA 1.1621, 0. DATA 0., 0.50088 DATA 0., 0.673 ' Fig. 5.8 DATA 0.01, 0.002, 0.002 ': mu, dx, dy DATA 8 ': NEP% DATA -1.14677, 0. ': dx, dy DATA -0.84808, 0. DATA 1.00417, 0. DATA -1.108, 0. DATA -1.040, 0. DATA -0.5597, 0. DATA -0.2894, 0. DATA 1.1127, 0. END FUNCTION dvdx (x, y, mu) ' ' evaluates partial derivative dV/dx in (x,y) ' dvdx = -(1 - mu) * (x - mu) / (r1(x, y, mu)) ^ 3 - mu *_ (x + 1 - mu) / (r2(x, y, mu) ^ 3) + x END FUNCTION FUNCTION dvdy (x, y, mu) ' ' evaluates partial derivative dV/dy in (x,y) ' dvdy = -(1 - mu) * y / (r1(x, y, mu) ^ 3)_ - mu * y / (r2(x, y, mu) ^ 3) + y END FUNCTION FUNCTION r1 (x, y, mu) ' ' computes distance from point (x,y) to first primary ' r1 = SQR((x - mu) ^ 2 + y ^ 2) END FUNCTION FUNCTION r2 (x, y, mu) ' ' compues distance from point (x,y) to second primary ' r2 = SQR((x + 1 - mu) ^ 2 + y ^ 2) END FUNCTION FUNCTION vxy (x, y, mu) ' ' evaluates potential energy function V(x,y) at point (x,y) ' vxy = (1 - mu) / r1(x, y, mu) + mu / r2(x, y, mu) _ + (x ^ 2 + y ^ 2) / 2 END FUNCTION SUB GRAPHSET(mu,xlps(),ylps()) xp2=mu-1 xp1=mu xL4=mu-0.5 yL4=0.8660254 SCREEN 18,2 __SCREENSET 1,0 WINDOW (-2.7,3.0)-(2.7,-0.5) VIEW (10,10)-(650,410) IX%=PMAP(xp1,0)/8+2 IY%=PMAP(-0.3,1)/16 LOCATE (IY%,IX%): PRINT "m1" CIRCLE(xp1,0),0.03,3 IX%=PMAP(xp2,0)/8+2 IY%=PMAP(-0.3,1)/16 LOCATE (IY%,IX%): PRINT "m2" CIRCLE(xp2,0),0.03,3 IX%=PMAP(xlps(1),0)/8+2 LOCATE (IY%,IX%): PRINT "L1" CIRCLE(xlps(1),0),0.01,3 IX%=PMAP(xlps(2),0)/8+2 LOCATE (IY%,IX%): PRINT "L2" CIRCLE(xlps(2),0),0.01,3 IX%=PMAP(xlps(3),0)/8+2 LOCATE (IY%,IX%): PRINT "L3" CIRCLE(xlps(3),0),0.01,3 IY%=PMAP(ylps(4),1)/16+1 IX%=PMAP(xlps(4),0)/8+3 LOCATE (IY%,IX%): PRINT "L4" CIRCLE(xlps(4),ylps(4)),0.01,3 ' Draw X and Y axes LINE (0.0,0)-(0.0,2.7),5 LINE (-2.7,0.0)-(2.7,0.0),5 LINE (xp1,0)-(xp1,0.1),5 LINE (xp2,0)-(xp2,0.1),5 LOCATE (26,44) END SUB SUB Nbranch(mu,x0,y0,k0,dx,dy,Nb%,xnxs(),ynxs()) r=sqr(dx^2+dy^2)*2 Ndiv%=1001 dth=3.141592653589793/(Ndiv%-1) Nb%=0 xpr=x0+dx ypr=y0 kpr=vxy(xpr,ypr,mu) thpr=0.0 IF(ABS(kpr-k0)<=1.E-5) THEN Nb%=1 xnxs(1)=xpr ynxs(1)=ypr ENDIF FOR I%=2 TO Ndiv% th=dth*(I%-1) x=r*cos(th)+x0 y=r*sin(th)+y0 k=Vxy(x,y,mu) IF (k-k0)*(kpr-k0) <=0. THEN thc=(th-thpr)*(k0-kpr)/(k-kpr)+thpr xc=r*cos(thc)+x0 yc=r*sin(thc)+y0 kc=Vxy(xc,yc,mu) Nb%=Nb%+1 xnxs(Nb%)=xc ynxs(Nb%)=yc err=ABS(1.0-k/k0)*100 ENDIF xpr=x ypr=y kpr=k thpr=th NEXT I% END SUB SUB LAGPS(mu,xlps(),ylps()) ' Find the locations of Lagrange points mum1=1-mu FOR I%=1 TO 3 IF I%=1 THEN x= 0.0 : h1=-1.0 : h2=1.0 ELSEIF I%=2 THEN x=-1.0 : h1=-1.0 : h2=-1.0 ELSE x= 1.0 : h1= 1.0 : h2=1.0 ENDIF c1=h1*mum1 c2=h2*mu FOR J%=1 TO 100 sa=x-mu sb=x+mum1 xnw=x-(x-c1/sa^2-c2/sb^2)/(2*(c1/sa^3+c2/sb^3)+1.0) IF(ABS(1.0-xnw/x)< 1.E-6) GOTO 200 x=xnw NEXT J% PRINT "Warning: not enough number of iterations!!!" 200: xlps(I%)=xnw ylps(I%)=0.0 NEXT I% xlps(4)=mu-0.5 ylps(4)=0.8660254 PRINT "No. of iterations=";J% END SUB