|
投稿者：SECOND
投稿日：2009年 6月 7日(日)00時37分5秒
|
|
|
|
! 立体曲線 f(x,y,z,t)の、変形指示MAT文 による描画。
!
! PLOT LINES :x1,y1,z1; x2,y2,z2 は、エラーですが、これを、picture 文に
! 置いて、draw picture with matrix したかの様に、描く。
!
! x1,y1,z1; x2,y2,z2 の両端 座標だけを、picture 文に渡して、
! draw … with … を実行し、座標のみの、変形、回転をする。
! x1,y1,z1; x2,y2,z2 のその出力から、x,y 成分だけの線を引く。
!（言い換えると、picture 文の中で、PLOT LINES をせず、出てから、PLOT する)
!-------------------------------
! ローレンツ方程式、３次元グラフ
!「例」描くグラフは、対流現象の近似式として知られるカオスです。
!-------------------------------
LET s=10
LET b=8/3
LET r=28
!-----ローレンツ方程式(カオス)
DEF dxdt(x,y )= -s*x +s*y ! (dx/dt)= -s*x +s*y
DEF dydt(x,y,z)= -x*z +r*x -y ! (dy/dt)= -x*z +r*x -y
DEF dzdt(x,y,z)= x*y -b*z ! (dz/dt)= x*y -b*z
SUB RungeKutta
LET kx1=dxdt(x,y )
LET ky1=dydt(x,y,z)
LET kz1=dzdt(x,y,z)
!
LET kx2=dxdt(x+kx1*dt/2, y+ky1*dt/2 )
LET ky2=dydt(x+kx1*dt/2, y+ky1*dt/2 ,z+kz1*dt/2)
LET kz2=dzdt(x+kx1*dt/2, y+ky1*dt/2 ,z+kz1*dt/2)
!
LET kx3=dxdt(x+kx2*dt/2, y+ky2*dt/2 )
LET ky3=dydt(x+kx2*dt/2, y+ky2*dt/2 ,z+kz2*dt/2)
LET kz3=dzdt(x+kx2*dt/2, y+ky2*dt/2 ,z+kz2*dt/2)
!
LET kx4=dxdt(x+kx3*dt, y+ky3*dt )
LET ky4=dydt(x+kx3*dt, y+ky3*dt ,z+kz3*dt)
LET kz4=dzdt(x+kx3*dt, y+ky3*dt ,z+kz3*dt)
!
LET x=x+(kx1+2*kx2+2*kx3+kx4)*dt/6
LET y=y+(ky1+2*ky2+2*ky3+ky4)*dt/6
LET z=z+(kz1+2*kz2+2*kz3+kz4)*dt/6
END SUB
!-----run
OPTION ANGLE DEGREES
DIM pV1(4), pV2(4), P3D(4,4), rotx(4,4)
MAT rotx=IDN
!
LET ax=-75 ! 原点を通り、画面の水平軸での回転（ x 軸とは、限らず。）
!
! 1, 0, 0, 0
! 0, cos(ax), 0, 0
! 0,-sin(ax), 1, 0
! 0, 0, 0, 1
!
LET rotx(2,2)=COS(ax)
LET rotx(3,2)=-SIN(ax)
!
LET xm=0
LET ym=20
LET h=40
SET WINDOW xm-h,xm+h,ym-h,ym+h
LET dt=.005 !sec. pitch time
!
FOR az=-45 TO 315 STEP 10 ! ｚ軸での回転。360 度、一回り。
!----
LET x=2
LET y=1
LET z=24
!----
LET t=0
IF -45< az THEN SET DRAW mode hidden
CLEAR
PLOT TEXT,AT xm+h*.2,ym+h*.9 :"lorenz ローレンツ方程式"
DO
IF 0< t THEN
!---3D 曲線(x,y,z)
DRAW line3D( bakx,baky,bakz, x,y,z) WITH ROTATE(az)*rotx
PLOT LINES: pV1(1),pV1(2); pV2(1),pV2(2)
ELSE
! ---X 軸
DRAW line3D( -30,0,0, 30,0,0) WITH ROTATE(az)*rotx
PLOT LINES: pV1(1),pV1(2); pV2(1),pV2(2)
PLOT TEXT,AT pV2(1),pV2(2) :"(X)"
!---Y 軸
DRAW line3D( 0,-30,0, 0,30,0) WITH ROTATE(az)*rotx
PLOT LINES: pV1(1),pV1(2); pV2(1),pV2(2)
PLOT TEXT,AT pV2(1),pV2(2) :"(Y)"
!---Z 軸
DRAW line3D( 0,0,-10, 0,0,50) WITH ROTATE(az)*rotx
PLOT LINES: pV1(1),pV1(2); pV2(1),pV2(2)
PLOT TEXT,AT pV2(1),pV2(2) :"(Z)"
END IF
!----
LET bakx=x
LET baky=y
LET bakz=z
CALL RungeKutta
LET t=t+dt
LOOP UNTIL 30< t ! 100 くらいが良いのかも… 遅くなる。
SET DRAW mode explicit
NEXT az
PICTURE line3D(x1,y1,z1, x2,y2,z2)
MAT P3D=TRANSFORM !← draw ･･･ with matrix で与えられた matrix
LET pV1(1)=x1
LET pV1(2)=y1
LET pV1(3)=z1 !入力 z1 座標は、出力 x1,y1 に反映、描画される。
LET pV1(4)=1
MAT pV1=pV1*P3D !出力 z1 座標pV1(3)は、描画しない。
LET pV2(1)=x2
LET pV2(2)=y2
LET pV2(3)=z2 !入力 z2 座標は、出力 x2,y2 に反映、描画される。
LET pV2(4)=1
MAT pV2=pV2*P3D !出力 z2 座標pV2(3)は、描画しない。
END PICTURE
END
|
|
|
|
投稿者：SECOND
投稿日：2009年 6月 9日(火)05時14分11秒
|
|
|
|
> No.407[元記事へ]
! 立体曲線 f(x,y,z,t)を、変形指示MAT文で、描く（２）
!
! 前回、くどい事をしていたようで、整理した。
! 配列ベクトル x,y,z 座標を入力に、単に、変形指示MAT文だけで変形回転する。
! x,y,z のその出力から、x,y 成分だけの PLOT。
!-------------------------------
! レスラー方程式 ３次元グラフ
! 非線形項が、z*x １つだけのカオス。形状:メビウスの帯
!-------------------------------
LET a=0.398
LET b=2
LET c=4
!-----レスラー方程式
SUB Dxyz( kx,ky,kz, x,y,z)
LET kx=-y-z ! (dx/dt)= -y -z
LET ky= x+a*y ! (dy/dt)= x +a*y
LET kz= b+z*(x-c) ! (dz/dt)= b +z*(x-c)
END SUB
SUB RungeKutta
CALL Dxyz( kx1,ky1,kz1, x,y,z)
CALL Dxyz( kx2,ky2,kz2, x+kx1*dt/2,y+ky1*dt/2,z+kz1*dt/2)
CALL Dxyz( kx3,ky3,kz3, x+kx2*dt/2,y+ky2*dt/2,z+kz2*dt/2)
CALL Dxyz( kx4,ky4,kz4, x+kx3*dt ,y+ky3*dt ,z+kz3*dt )
LET x=x+(kx1+2*kx2+2*kx3+kx4)*dt/6
LET y=y+(ky1+2*ky2+2*ky3+ky4)*dt/6
LET z=z+(kz1+2*kz2+2*kz3+kz4)*dt/6
END SUB
!-----run
OPTION CHARACTER byte !ウムラウト"o" =CHR$(246) をエラーにさせない。
OPTION ANGLE DEGREES
DIM pV(4), P3D(4,4), rotx(4,4)
!
MAT rotx=IDN ! 単位行列
LET ax=-55 ! 原点を通り、画面の水平軸での回転（ x 軸とは、限らず。）
!
!(x,y,z,1)| 1, 0, 0, 0 | PLOT 座標ベクトルは、行ベクトル。
! | 0, cos(ax), 0, 0 | ４列目の１は、拡大率の逆数であるが、
! | 0,-sin(ax), 1, 0 | draw 文で、描画しないので、無視して可。
! | 0, 0, 0, 1 |
!
LET rotx(2,2)=COS(ax)
LET rotx(3,2)=-SIN(ax)
!
LET xm=0
LET ym=2
LET h=7
SET WINDOW xm-h,xm+h,ym-h,ym+h
LET dt=.05 ! sec. pitch time
LET SS=35 ! ｚ軸回り、開始角度
LET EE=SS-360 ! +360：左回転 -360：右回転
!
FOR az=SS TO EE STEP SGN(EE-SS)*10 ! ｚ軸、一回り360 度
MAT P3D=ROTATE(az)*rotx
!----
IF az<>SS THEN SET DRAW mode hidden
CLEAR
SET TEXT font "Courier",11
PLOT TEXT,AT xm+h*.2,ym+h*.9 :"R"& CHR$(246)& "ssler"
SET TEXT font "標準ゴシック",11
PLOT TEXT,AT xm+h*.45,ym+h*.9 :"レスラー方程式" ! PEN-off
!---座標軸
CALL axes3D( -6,0,0, 6,0,0, "(X)" )
CALL axes3D( 0,-6,0, 0,6,0, "(Y)" )
CALL axes3D( 0,0,-2, 0,0,8, "(Z)" )
!---3D 曲線
LET x=-3
LET y=0
LET z=0
FOR t=0 TO 300 STEP dt
CALL line3D(x,y,z)
CALL RungeKutta
NEXT t
SET DRAW mode explicit
NEXT az
SUB axes3D(x1,y1,z1, x2,y2,z2, a$ )
CALL line3D(x1,y1,z1)
CALL line3D(x2,y2,z2)
PLOT TEXT,AT pV(1),pV(2) :a$ ! PEN-off
END SUB
SUB line3D(x,y,z)
LET pV(1)=x
LET pV(2)=y
LET pV(3)=z !入力 z 座標は、出力 x,y に反映、描画される。
MAT pV=pV*P3D
PLOT LINES: pV(1),pV(2); !出力 z 座標pV(3)は、不要。 PEN-on
END SUB
END
|
|
|
|
投稿者：SECOND
投稿日：2009年 6月12日(金)05時12分43秒
|
|
|
|
> No.407[元記事へ]
!先の投稿、ローレンツ方程式で、DEF dxdt()、DEF・・・の引数に脱落が
!ありました。すみません。 （ 現在、訂正済み。）
!-------------------------------
! 立体曲線 f(x,y,z,t)を、変形指示MAT文で、描く（３）
!-------------------------------
OPTION ARITHMETIC NATIVE
OPTION ANGLE DEGREES
DIM pV(4), P3D(4,4), rotx(4,4), LH(6), copy(0 TO 100000, 3)
MAT rotx=IDN ! 単位行列
!
!-----
! 電子回路で 作られたカオス。(添付図参照)
LET t$="Double Scroll Attractor ダブル スクロール アトラクタ"
! 変曲点(±Bp)のある 負性抵抗回路への 電流＝ih( 電圧 )
DEF ih(x)= g0*x+(g1-g0)*( ABS(x+Bp)-ABS(x-Bp) )/2
SUB Dxyz( kx,ky,kz, x,y,z)
LET kx= (gR*(y-x) -ih(x))/C1 ! (d vC1/dt)= (gR*(vC2-vC1)-ih(vC1))/C1
LET ky= (gR*(x-y) +z )/C2 ! (d vC2/dt)= (gR*(vC1-vC2)+iL )/C2
LET kz= -y/L ! (d iL /dt)= -vC2/L
END SUB
LET C1= 1/9 !コンデンサー （パラメーター C1~Bp）
LET C2= 1 !コンデンサー
LET L= 1/7 !インダクター
LET gR= 0.7 !Ｒのコンダクタンス（アドミタンスの実数部）
LET g0=-0.5 !負性抵抗回路の微分コンダクタンス（～＜-Bp +Bp＜～)
LET g1=-0.8 !負性抵抗回路の微分コンダクタンス（ -Bp＜～＜+Bp )
LET Bp= 1 !負性抵抗回路の、変曲点の電圧 絶対値。
!
LET x=-1e-7 !初期値 x,y,z
LET y= 1e-7
LET z= 1e-7
DATA -3,3, -3,3, -2.5,3 !座標軸の両端 xL,xH, yL,yH, zL,zH
MAT READ LH
LET xm=0 !画面中心 xm,ym
LET ym=.5 !
LET h=3.5 !画面片幅 ±h
LET dt=.01 ! pitch time
LET t99=120 ! close time
LET ax=-85 !ｚ軸を、画面の水平軸で回転、傾ける角度
LET SS=-25 !ｚ軸での、回転 開始角度
LET EE=SS+360 ! 回転 終了角度 +360：左回転 -360：右回転
CALL graph3D
!-----
SUB RungeKutta
CALL Dxyz( kx1,ky1,kz1, x,y,z)
CALL Dxyz( kx2,ky2,kz2, x+kx1*dt/2,y+ky1*dt/2,z+kz1*dt/2)
CALL Dxyz( kx3,ky3,kz3, x+kx2*dt/2,y+ky2*dt/2,z+kz2*dt/2)
CALL Dxyz( kx4,ky4,kz4, x+kx3*dt ,y+ky3*dt ,z+kz3*dt )
LET x=x+(kx1+2*kx2+2*kx3+kx4)*dt/6
LET y=y+(ky1+2*ky2+2*ky3+ky4)*dt/6
LET z=z+(kz1+2*kz2+2*kz3+kz4)*dt/6
END SUB
SUB graph3D
!-----run
!ｚ軸を、画面の水平軸で回転、傾ける行列 rotx
!(x,y,z,1)| 1, 0, 0, 0 | 十進の、PLOTベクトルは、行ベクトル。
! | 0, cos(ax), 0, 0 | ４列目の１は、拡大率の逆数であるが、
! | 0,-sin(ax), 1, 0 | draw~with~ で効果する。
! | 0, 0, 0, 1 | このプログラムでは、無視して可。
LET rotx(2,2)=COS(ax)
LET rotx(3,2)=-SIN(ax)
!
SET WINDOW xm-h,xm+h,ym-h,ym+h
!
FOR az=SS TO EE STEP SGN(EE-SS)*10 ! ｚ軸で、一回り360 度
MAT P3D=ROTATE(az)*rotx
!----
IF az<>SS THEN SET DRAW mode hidden
CLEAR
PLOT TEXT,AT xm-h*.8,ym+h*.9 :t$ ! PEN-off
!---座標軸
CALL axes3D( LH(1),0,0, LH(2),0,0, "(X)" )
CALL axes3D( 0,LH(3),0, 0,LH(4),0, "(Y)" )
CALL axes3D( 0,0,LH(5), 0,0,LH(6), "(Z)" )
!---3D 曲線
IF az=SS THEN
LET n=0
FOR t=0 TO t99 STEP dt
LET copy(n,1)=x
LET copy(n,2)=y ! １回目(開始角度)で、３Ｄ記録を撮る。
LET copy(n,3)=z
CALL line3D(x,y,z)
CALL RungeKutta
LET n=n+1
NEXT t
ELSE
FOR n=0 TO n-1
LET x=copy(n,1)
LET y=copy(n,2) ! ２回目以降は、記録の再生で、高速描画。
LET z=copy(n,3)
CALL line3D(x,y,z)
NEXT n
END IF
SET DRAW mode explicit
NEXT az
END SUB
SUB axes3D(x1,y1,z1, x2,y2,z2, a$ )
CALL line3D(x1,y1,z1)
CALL line3D(x2,y2,z2)
PLOT TEXT,AT pV(1),pV(2) :a$ ! PEN-off
END SUB
SUB line3D(x,y,z)
LET pV(1)=x
LET pV(2)=y
LET pV(3)=z !入力 z 座標は、出力 x,y に反映、描画される。
MAT pV=pV*P3D
PLOT LINES: pV(1),pV(2); !出力 z 座標pV(3)は、不要。 PEN-on
END SUB
END
|
|
|
戻る
