[VA] SRC #012f - Then and Now: Angle

[J-F Garnier]

(05-08-2023 08:40 AM)Fernando del Rey Wrote:  Thanks, J-F, your two improvements in post #4 make a significant difference in terms of execution speed.

I have tried to optimize speed by using an adaptative H for the Runge-Kutta, with larger H when the absolute value of the acceleration vector A is low, and smaller H when acceleration is high. But I did not get any maningful results out of this idea.

There is an easy way to shorten the search for the root, since we are seeking for only 5 significant digits, with this well known trick:

      145 IF ABS(L)<1.E-5 THEN L=0 ! Close enough from the target !

Try it !

(05-08-2023 10:08 PM)Valentin Albillo Wrote:  I'll wait some extra days before posting my original solution (concocted almost a year ago), just in case that someone else wants to post something, you optimize your solution even further and/or J-F posts his very own solution.

So it's time for me to release the latest optimizations I was able to imagine:
- check each trajectory for both the target (4,1) and the symmetric (4,-1) positions, this let us search for positive angles only,
- rewrite the FNA function without loop and without exponentiation (a pure programming optimization),
- replace INTEGER variables with REAL (no speed gain to use INTEGER in 71B BASIC).

The modified program becomes (added/modified lines in bold):

10 DESTROY ALL @ OPTION BASE 1 @ RADIANS
20 COMPLEX P,V,A,S(6),Z,M1,M2,M3,M4,K1,K2,K3,K4
30 REAL I,J,Z0,D,D5,D6,H,L,T,X,E(50),U(50)
40 DATA (2,-1),(2,0),(2,1),(3,-1),(3,0),(3,1),(4,1)
50 ! Function to calculate acceleration vector at point P
60 ! DEF FNA(P) @ A=0 @ FOR I=1 TO 6 @ A=A+(S(I)-P)/ABS(S(I)-P)^3 @ NEXT I @ FNA=A @ END DEF
61 DEF FNA(P) @ X=ABS(S(1)-P) @ A=(S(1)-P)/(X*X*X) @ X=ABS(S(2)-P) @ A=A+(S(2)-P)/(X*X*X)
62 X=ABS(S(3)-P) @ A=A+(S(3)-P)/(X*X*X) @ X=ABS(S(4)-P) @ A=A+(S(4)-P)/(X*X*X)
63 X=ABS(S(5)-P) @ A=A+(S(5)-P)/(X*X*X) @ X=ABS(S(6)-P) @ A=A+(S(6)-P)/(X*X*X)
64 FNA=A @ END DEF
70 ! Function to calculate trajectory using Runge-Kutta method
80 DEF FNP(D) @ T=0 @ V=(COS(D),SIN(D)) @ P=(0,0)
90 M1=H*V @ K1=H*FNA(P) @ M2=H*(V+K1/2) @ K2=H*FNA(P+M1/2)
100 M3=H*(V+K2/2) @ K3=H*FNA(P+M2/2) @ M4=H*(V+K3) @ K4=H*FNA(P+M3)
110 P=P+(M1+M2+M2+M3+M3+M4)/6 @ V=V+(K1+K2+K2+K3+K3+K4)/6
120 T=T+H @ IF REPT(P)<Z0 AND T<T9 THEN 90
130 IF T>=T9 THEN L=0 @ GOTO 150 ! If time is too long, exit
135 IF IMPT(P)>0 THEN D6=D @ L=IMPT(P)-IMPT(Z) ELSE D6=-D @ L=IMPT(P)+IMPT(Z)
140 L=L-(REPT(P)-Z0)*(IMPT(V)/REPT(V)) @ T=T-(REPT(P)-Z0)/REPT(V)
145 IF ABS(L)<1.E-5 THEN L=0 ! Close enough from the target !
150 FNP=L @ END DEF ! L=distance from target vertical line, T=travel time
160 ! Main program SRC #12f - with JFG's modif
170 T0=TIME @ T5=MAXREAL @ T9=2.5 ! T9 is max trajectory time
180 FOR I=1 TO 6 @ READ S(I) @ NEXT I @ READ Z @ Z0=REPT(Z) ! Read star and target positions
190 FOR D5=0 TO .75 STEP .0125 ! Angle search range and step
200 H=.025 @ DISP D5 ! Delta-T for coarse trajectory
210 IF J=0 THEN 230
220 IF D5<ABS(E(J))+.02 THEN 320 ! If too close from a root angle already found, skip
230 IF ABS(FNP(D5))>.5 THEN 320 ! If too far from target, skip
235 IF ABS(V)>2 THEN 320 ! If abnormal final speed, skip
240 IF T>=T9 THEN 320 ! If time is too long, skip
250 H=.01 @ DISP D6;"*" ! Delta-T for fine trajectory
260 D0=FNROOT(D6,D6,FNP(FVAR)) ! Fine search of root
270 IF T>=T9 THEN 320 ! If time is too long, skip
280 IF J=0 THEN 300
290 IF ABS(D0-E(J))<1.E-11 THEN 320 ! If root angle was already found, skip
300 J=J+1 @ E(J)=D0 @ U(J)=T ! Save root angle and trajectory time
310 IF T<T5 THEN T5=T @ J5=J ! Look for shortest trajectory time
320 NEXT D5 @ T0=TIME-T0
330 DISP "Angle= ";E(J5);"Time= ";T5;T0

The result is:

Angle= -.274882047936 Time= 1.91608959636 (96.61s)

obtained in about 97s on my Emu71/Win system with ~1400 speed factor,
i.e. about 38 h on a real HP-71B. Not yet a run overnight, but very possible over a week-end.

The 6 trajectories:
    Angle             Time
-.196990679676 1.99998010763
-.274882047936 1.91608959636 
 .328560456480 2.0237867079
-.461225625588 2.1526452933
 .551285709122 1.91762496347
 .732203223220 2.06532676912

J-F
Edit: typos.