|
|
! 立体曲線 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)
LET pV(4)=1 ! shift()*rotate()* … で必要。
MAT rotx=IDN
!
!-----
LET t$="ヤリイカの巨大神経（カオス） 数学モデル"
LET t2$="Hodgkin-Huxley ホジキン－ハクスレイ方程式 のストレンジ アトラクタ"
DEF I(t)= I0+A*SIN( 360*f*t ) ! 膜電流 +外部入力
LET t3$="膜電位 V →( X)"
! (d V/dt)= I(t)-120*m^3*h*(V-115) -40*n^4*(V+12) -.24*(V-10.613)
!
LET t4$="ナトリウム 活性化変数( 0< m< 1) →( Y)"
! (d m/dt)= .1*(25-V)/(EXP((25-V)/10)-1)*(1-m) -4*EXP(-V/18)*m
!
LET t5$="ナトリウム不活性化変数( 0< h< 1) →( Z)"
! (d h/dt)= .07*EXP(-V/20)*(1-h) -1/(EXP((30-V)/10)+1)*h
! カリウム活性化変数( 0< n< 1) →描画しない。
! (d n/dt)= .01*(10-V)/(EXP((10-V)/10)-1)*(1-n) -.125*EXP(-V/80)*n
SUB Dxyzn( kx,ky,kz,kn, x,y,z,n)
LET kx= I(t)-120*y^3*z*(x-115) -40*n^4*(x+12) -.24*(x-10.613) ! (d V/dt)
LET ky= .1*(25-x)/(EXP((25-x)/10)-1)*(1-y) -4*EXP(-x/18)*y ! (d m/dt)
LET kz= .07*EXP(-x/20)*(1-z) -1/(EXP((30-x)/10)+1)*z ! (d h/dt)
LET kn= .01*(10-x)/(EXP((10-x)/10)-1)*(1-n) -.125*EXP(-x/80)*n! (d n/dt)
END SUB
LET I0=20 !膜電流 パラメーター
LET A=40
LET f=.3000001
!
LET x=6.24 !初期値 x,y,z,n
LET y=.0761
LET z=.301
LET n=.519
!
DATA -15,100, -.3,1, -.04,.49 !座標軸の両端座標 xL,xH, yL,yH, zL,zH
MAT READ LH
LET Sx=1 !スケール倍率 Sx,Sy,Sz
LET Sy=50
LET Sz=200
!
LET xm=40 !画面中心 xm,ym
LET ym=52
LET hw=80 !画面幅/2 ±hw
LET dt=.05 ! pitch time
LET t99=500 ! close time
LET zox=40 !ｚ軸と、平行な回転軸へのオフセットｘｙ
LET zoy=20
LET ax=-75 !回転軸を、画面の水平軸で倒し、傾ける角度
LET SS=-25 !回転 開始角度
LET EE=SS+360 !回転 終了角度 +360：左回転 -360：右回転
CALL graph3D
!-----
SUB RungeKutta
CALL Dxyzn( kx1,ky1,kz1,kn1, x,y,z,n)
CALL Dxyzn( kx2,ky2,kz2,kn2, x+kx1*dt/2,y+ky1*dt/2,z+kz1*dt/2,n+kn1*dt/2)
CALL Dxyzn( kx3,ky3,kz3,kn3, x+kx2*dt/2,y+ky2*dt/2,z+kz2*dt/2,n+kn2*dt/2)
CALL Dxyzn( kx4,ky4,kz4,kn4, x+kx3*dt ,y+ky3*dt ,z+kz3*dt ,n+kn3*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
LET n=n+(kn1+2*kn2+2*kn3+kn4)*dt/6
END SUB
SUB graph3D
!-----run
!回転軸を、画面の水平軸で倒し、傾ける行列 rotx
!(x,y,z,1)| 1, 0, 0, 0 | 十進の、PLOTベクトルは、行ベクトル。
! | 0, cos(ax), 0, 0 | ４列目の１は、拡大率の逆数、draw文 で効果。
! | 0,-sin(ax), 1, 0 | 直接の plot lines では、1 にリセットされる。
! | 0, 0, 0, 1 | 行列の計算では、平行移動を、可能にする。
LET rotx(2,2)=COS(ax) !このプログラムでは、SHIFT(,) で必要。
LET rotx(3,2)=-SIN(ax)
!
SET WINDOW xm-hw,xm+hw,ym-hw,ym+hw
!
FOR az=SS TO EE STEP SGN(EE-SS)*10 ! 回転軸で、10 度づつ回す。
MAT P3D=SHIFT(-zox,-zoy)*ROTATE(az)*SHIFT(zox,zoy)*rotx
!----
IF az<>SS THEN SET DRAW mode hidden
CLEAR
PLOT TEXT,AT xm-hw*.9,ym+hw*.90 :t$
PLOT TEXT,AT xm-hw*.9,ym+hw*.83 :t2$
PLOT TEXT,AT xm-hw*.9,ym-hw*.92 :t5$
PLOT TEXT,AT xm-hw*.9,ym-hw*.99 :t4$& " "& t3$ ! PEN-off
!---座標軸
CALL axes3D( LH(1),0,0, LH(2),0,0, STR$(LH(2))& "( X)" )
CALL axes3D( 0,LH(3),0, 0,LH(4),0, STR$(LH(4))& "( Y)" )
CALL axes3D( 0,0,LH(5), 0,0,LH(6), STR$(LH(6))& "( Z)" )
!---3D 曲線
IF az=SS THEN
LET ci=0
FOR t=0 TO t99 STEP dt
LET copy(ci,1)=x
LET copy(ci,2)=y ! １回目(開始角度)で、３Ｄ記録を撮る。
LET copy(ci,3)=z
! PRINT x;y;z;n ! データーを保存したい時。
CALL line3D(x,y,z)
CALL RungeKutta
LET ci=ci+1
NEXT t
ELSE
FOR ci=0 TO ci-1 ! ２回目以降は、記録の再生で、高速描画。
CALL line3D( copy(ci,1),copy(ci,2),copy(ci,3) )
NEXT ci
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*Sx !描画目盛は、全方向等しくないと、回転で、形が保てない。
LET pV(2)=y*Sy !スケール Sx,Sy,Sz の違いは、入力の倍率として、行なう。
LET pV(3)=z*Sz !入力 z 座標は、出力 x,y に反映、描画される。
MAT pV=pV*P3D
PLOT LINES: pV(1),pV(2); ! PEN-on
END SUB
END
|
|