10 '-----ピストン加速度の近似_サブ------
20 CLS 3:SCREEN 3
30 X=100:Y=200 : TIME=1000
40 PAI=3.141592653#
50 RAD=PAI/180
60 ZOOM=80
70 A=.8#
80 LINE(X-50,Y)-(X+400,Y)
90 LINE(X,Y-150)-(X,Y+150)
100 '-----厳密解-----
110 FOR T=0 TO 360 STEP 1
120 THETA=T*RAD
130 V1=-(COS(THETA)+(((A*COS(2*THETA))+(A^3*SIN(THETA)^4))/(SQR(1-(A*SIN(THETA))^2)^3)))
140 IF T=0 THEN 150 ELSE 160
150 PSET(X+T,Y-V1*ZOOM):GOTO 170
160 LINE -(X+T,Y-V1*ZOOM)
170 FOR W=0 TO TIME STEP 1 : NEXT W
180 NEXT T
190 '------近似解------
200 B=0
210 B(0)=COS(THETA)
220 B(1)=(4*(((1/4)*A)+((1/16)*A^3)+((15/512)*A^5)+((35/2048)*A^7)+((735/65536)*A^9)+((2079/262144)*A^11)+((99099/16777216)*A^14)))
230 B(2)=(16*(((-1/64)*A^3)-((3/256)*A^5)-((35/4096)*A^7)-((105/16384)*A^9)-((10395/2097152)*A^11)-((33033/8388608)*A^14)))
240 B(3)=(36*(((1/512)*A^5)+((5/2048)*A^7)-((315/131072)*A^9)-((1155/524288)*A^11)+((33033/16777216)*A^14)))
250 B(4)=(64*(((-5/16384)*A^7)-((35/65536)*A^9)-((693/1048576)*A^11)-((3003/4194304)*A^14)))
260 B(5)=(100*(((7/131072)*A^9)-((63/524288)*A^11)+((3003/16777216)*A^14)))
270 B(6)=(144*(((-21/2097152)*A^11)-((231/8388608)*A^14)))
280 B(7)=(196*(((-33/16777216)*A^14)))
290 N=0
300 LOCATE 30,2 : INPUT "nの値を入力せよ";NN
310 FOR II=0 TO NN STEP 1
320 FOR T=0 TO 359 STEP 1
330 THETA1=T*RAD
340 THETA2=(T+1)*RAD
350 IF II=0 THEN 
360 III=0.5
370 Y1=-B(II)*COS(2*III*THETA1)*ZOOM
380 Y2=-B(II)*COS(2*III*THETA2)*ZOOM
390 LINE (X+T,Y-Y1)-(X+T+1,Y-Y2),4
400 ELSE
410 III=II
420 Y1=-B(II)*COS(2*III*THETA1)*ZOOM
430 Y2=-B(II)*COS(2*III*THETA2)*ZOOM
440 LINE (X+T,Y-Y1)-(X+T+1,Y-Y2),4
450 END IF
460 NEXT T
470 LOCATE 15,3 : PRINT "Enter キーを押してください"
480 IF INKEY$ = "" THEN
490 GOTO 480
500 ELSE
510 GOTO 520 							
520 END IF
530 NEXT II
540 FOR T=0 TO 360
550 C=0
560 FOR I=1 TO NN
570 THETA=T*RAD
580 IF I=1 THEN 
590 C=-COS(THETA)-(B(1)*COS(2*I*THETA)) 
600 ELSE 
610 C=C-(B(I)*COS(2*I*THETA))
620 END IF
630 NEXT I
640 IF T=0 THEN 
650 PSET(X+T,Y-C*ZOOM),4
660 ELSE
670 LINE-(X+T,Y-C*ZOOM),3
680 END IF
690 NEXT T
700 END