teacup. [ 掲示板 ] [ 掲示板作成 ] [ 有料掲示板 ] [ ブログ ]

 投稿者
  題名
  内容 入力補助画像・ファイル<IMG>タグが利用可能です。(詳細)
    
 URL
[ ケータイで使う ] [ BBSティッカー ] [ 書込み通知 ] [ 検索 ]

スレッド一覧

  1. スレッドが使えます(2)
  2. Paract BASIC(22)
  3. Amusement_Program(10)
  4. 改修予定のJIS非互換(3)
  5. 複数ページ長編プログラム(新規投稿)(16)
  6. 「十進BASIC第2掲示板」投稿記事リスト(17)
  7. Full BASIC互換ライブラリ(8)
  8. 「十進BASIC掲示板過去ログ」インデックス(トピック)(17)
  9. 人の色覚の数理(14)
  10. 「十進BASIC掲示板過去ログ」インデックス(ツリー)(91)
スレッド一覧(全10)  他のスレッドを探す  スレッド作成

*掲示板をお持ちでない方へ、まずは掲示板を作成しましょう。無料掲示板作成


Re: Yahoo!ジオシティーズ サービス終了について

 投稿者:白石 和夫  投稿日:2018年10月 9日(火)14時52分2秒
返信・引用
  > No.4565[元記事へ]

忍者ホームページに移設しました。
http://decimalbasic.ninja-web.net/
不具合など見つけたらお知らせください。
 
 

Yahoo!ジオシティーズ サービス終了について

 投稿者:marimo  投稿日:2018年10月 7日(日)18時19分29秒
返信・引用
  Yahoo!ジオシティーズが2019年3月末をもってサービス終了という発表がありました 。
<参考>
サービス終了のお知らせ - Yahoo!ジオシティーズ https://info-geocities.yahoo.co.jp/close/index.html

Mathematics Education http://www.geocities.jp/thinking_math_education/
で公開されてるページが2019年3月末でアクセスできなくなるかと思われます。
今後他のサービスなどへの移行は行われるのでしょうか?

個人的にWebヘルプを良く使ってるので気になります。
(この話、既出だったらすいません)
 

SAMPLE\Collatz.basにバグがあります。

 投稿者:hayashi  投稿日:2018年 9月 4日(火)15時22分5秒
返信・引用
  SAMPLE\Collatz.basにバグがあります。
---------------------------------
REM コラッツ予想
REM 任意の自然数に対し,
REM 偶数のとき,2で割り,
REM 奇数のとき,3倍して1を加える
REM という操作を繰り返すと,
REM 有限回で1に到達すると予想されている。
REM 証明はまだない。
OPTION ARITHMETIC RATIONAL
INPUT n
DO
   IF MOD(n,2)=0 THEN
      LET n=n/2
   ELSE
      LET n=3*n+1
      PRINT n
   END IF
LOOP UNTIL n=1
END

-------------------
      PRINT n
   END IF
この2行を入れ替えないと
1ではなく16で止まってしまいます。
 

リンク切れのデータ

 投稿者:田中  投稿日:2018年 8月11日(土)16時42分38秒
返信・引用
  時々、拝見しております。もうすでに回答されていることかと思いますが、
 山中和義 氏の
  センター数学(BASIC)による「ゲームプログラミングの宝箱」
2016年頃に、上記の書庫がZIPで、ダウンロード可能だったと思います。当方もその項目リストだけは、プリントしたのですが、本体のプログラムは、紛失してしまったか、不明になってしまいました。
 これの所在は、どこかにないでしょうか。ときどきリンクして紹介されていますが、クリックしてもリンク切れになっていますので、研究用に参考にしたく思っています。よろしくお願いします。
 

Re: 座標軸の太さ

 投稿者:白石 和夫  投稿日:2018年 6月 2日(土)17時46分53秒
返信・引用
  > No.4561[元記事へ]

nagramさんへのお返事です。

独自拡張なのでこれが仕様と考えてください。
独自の軸描画が必要であれば、Libraryフォルダのaxes.lib, grid.libに手をいれて使ってください。


> 十進BASICの独自拡張である軸・格子を描く絵定義 axes, grid が、set line width による設定の影響を受けます。
> これは仕様でしょうか、バグでしょうか。
>
> SET WINDOW -5,5,-5,5
> SET LINE width 3
> DRAW axes
> pause
> CLEAR
> SET LINE width 5
> DRAW grid
> END
 

座標軸の太さ

 投稿者:nagram  投稿日:2018年 6月 2日(土)14時15分53秒
返信・引用
  十進BASICの独自拡張である軸・格子を描く絵定義 axes, grid が、set line width による設定の影響を受けます。
これは仕様でしょうか、バグでしょうか。

SET WINDOW -5,5,-5,5
SET LINE width 3
DRAW axes
pause
CLEAR
SET LINE width 5
DRAW grid
END
 

Re: オーディオスペクトラム画像

 投稿者:しばっち  投稿日:2018年 5月20日(日)20時15分10秒
返信・引用
  > No.4559[元記事へ]

続き

SUB POWERDISK(RR(),XS,XE,YS,YE,COL,MAXCOL,LEV,MODE$)
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET P=10*LOG10(2^(SAMPLEBIT-1))
   LET YS=-P
   LET YE=P
   LET XE=P
   LET XS=-P
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN) !'ハニング窓関数を掛ける
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI) !'FFT
   LET T=NN/LEV/2
   FOR I=0 TO NN/2-1 STEP T
      LET R=0
      FOR P=0 TO T-1
         LET R=R+10*LOG10(FR(I+P)^2+FI(I+P)^2+1) !'平均値
      NEXT P
      LET R=R/2/T
      LET TH1=I/(NN/2-1)*360
      LET TH2=TH1+(T-NN/2/32)/(NN/2)*360
      SELECT CASE MODE$
      CASE "MODE1"
         FOR Y=R TO YE/16  STEP -YE/32
            IF Y<MAX(YS,YE)/2 THEN CALL DISK(0,0,Y+YE/32,COL,TH1,TH2) ELSE CALL DISK(0,0,Y+YE/32,MAXCOL,TH1,TH2)
         NEXT Y
         CALL DISK(0,0,YE/16,7,0,360)
      CASE "MODE2"
         IF MOD(R,YE/8)<YE/16 THEN LET X=0 ELSE LET X=1
         FOR Y=INT(R/(YE/16))*(YE/16) TO YE/16 STEP -YE/16
            LET X=1-X
            IF X=1 THEN
               IF Y<MAX(YS,YE)/2 THEN CALL DISK(0,0,Y+YE/16,COL,TH1,TH2) ELSE CALL DISK(0,0,Y+YE/16,MAXCOL,TH1,TH2)
            ELSE
               CALL DISK(0,0,Y+YE/16,0,TH1,TH2)
            END IF
         NEXT Y
         CALL DISK(0,0,YE/16,7,0,360)
      END SELECT
   NEXT I
END SUB

SUB VOLUMEBAR(RR(),XS,XE,YS,YE,COL,MAXCOL,MODE$,DIRECTION$) !'音量バー
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET P=10*LOG10(2^(SAMPLEBIT-1))
   SELECT CASE DIRECTION$
   CASE "BOTTOM"
      LET YS=0
      LET YE=P
      LET XE=1
      LET XS=0
   CASE "TOP"
      LET YS=P
      LET YE=0
      LET XS=0
      LET XE=1
   CASE "LEFT"
      LET XS=0
      LET XE=P
      LET YS=0
      LET YE=1
   CASE "RIGHT"
      LET XS=P
      LET XE=0
      LET YS=0
      LET YE=1
   END SELECT
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN)
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
   LET R=0
   FOR I=0 TO NN/2-1
      LET R=R+10*LOG10(FR(I)^2+FI(I)^2+1) !'平均値
   NEXT I
   LET R=R/(NN/2)
   SELECT CASE DIRECTION$
   CASE "LEFT","RIGHT"
      SELECT CASE MODE$
      CASE "MODE1"
         IF R>MAX(XS,XE)/32 THEN
            FOR X=0 TO R STEP MAX(XS,XE)/32
               IF X<MAX(XS,XE)*.8 THEN CALL BOXFULL(X,0,X+MAX(XS,XE)/32,1,COL) ELSE CALL BOXFULL(X,0,X+MAX(XS,XE)/32,1,MAXCOL)
            NEXT X
         END IF
      CASE "MODE2"
         IF R>MAX(YS,YE)/16 THEN
            FOR X=0 TO R STEP MAX(XS,XE)/16
               IF X<MAX(XS,XE)*.8 THEN CALL BOXFULL(X,0,X+MAX(XS,XE)/16*.8,1,COL) ELSE CALL BOXFULL(X,0,X+MAX(XS,XE)/16*.8,1,MAXCOL)
            NEXT X
         END IF
      END SELECT
   CASE "TOP","BOTTOM"
      SELECT CASE MODE$
      CASE "MODE1"
         IF R>MAX(YS,YE)/32 THEN
            FOR Y=0 TO R STEP MAX(YS,YE)/32
               IF Y<MAX(YS,YE)*.8 THEN CALL BOXFULL(0,Y,1,Y+MAX(YS,YE)/32,COL) ELSE CALL BOXFULL(0,Y,1,Y+MAX(YS,YE)/32,MAXCOL)
            NEXT Y
         END IF
      CASE "MODE2"
         IF R>MAX(YS,YE)/16 THEN
            FOR Y=0 TO R STEP MAX(YS,YE)/16
               IF Y<MAX(YS,YE)*.8 THEN CALL BOXFULL(0,Y,1,Y+MAX(YS,YE)/16*.8,COL) ELSE CALL BOXFULL(0,Y,1,Y+MAX(YS,YE)/16*.8,MAXCOL)
            NEXT  Y
         END IF
      END SELECT
   END SELECT
END SUB

SUB VOLUMECIRCLE(RR(),XS,XE,YS,YE,COL,MAXCOL) !'音量円形
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET YS=-10*LOG10(2^(SAMPLEBIT-1))
   LET YE=-YS
   LET XS=-10*LOG10(2^(SAMPLEBIT-1))
   LET XE=-XS
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN)
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
   LET R=0
   FOR I=0 TO NN/2-1
      LET R=R+10*LOG10(FR(I)^2+FI(I)^2+1) !'平均値
   NEXT I
   LET R=R/(NN/2)/2
   IF COUNT=1 OR MOD(COUNT,INT(FRAMERATE/2))=0 THEN LET TMAX=0
   LET TMAX=MAX(TMAX,R)
   CALL CIRCLE(0,0,TMAX,MAXCOL)
   IF MOD(R,XE/8)<XE/16 THEN LET Y=0 ELSE LET Y=1
   FOR X=INT(R/(XE/16))*(XE/16) TO XE/16 STEP -XE/16
      LET Y=1-Y
      IF Y=1 THEN
         CALL CIRCLEFULL(0,0,X,COL)
      ELSE
         CALL CIRCLEFULL(0,0,X,0)
      END IF
   NEXT X
   CALL CIRCLEFULL(0,0,XE/16,7)
END SUB

SUB XYPLOT(RR(),XS,XE,YS,YE,COL) !'スペクトルと位相
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET YS=-PI
   LET YE=PI
   LET XS=0
   LET XE=10*LOG10(2^(SAMPLEBIT-1))
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN)
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
   SET COLOR COL
   FOR I=0 TO NN/2-1
      IF FR(I)<>0 OR FI(I)<>0 THEN PLOT POINTS:10*LOG10(FR(I)^2+FI(I)^2+1),ANGLE(FR(I),FI(I))
   NEXT I
END SUB
END

EXTERNAL SUB BOXFULL(X1,Y1,X2,Y2,C)
OPTION ARITHMETIC NATIVE
SET AREA COLOR C
PLOT AREA:X1,Y1;X2,Y1;X2,Y2;X1,Y2
END SUB

EXTERNAL SUB CIRCLEFULL(X,Y,RR,C)
SET COLOR C
DRAW DISK WITH SCALE(RR)*SHIFT(X,Y)
END SUB

EXTERNAL SUB CIRCLE(X,Y,R,C)
SET COLOR C
DRAW CIRCLE WITH SCALE(R)*SHIFT(X,Y)
END SUB

EXTERNAL SUB DISK(XX,YY,R,C,S,T)
OPTION ANGLE DEGREES
OPTION BASE 0
DIM X(361),Y(361)
MAT X=(XX)*CON
MAT Y=(YY)*CON
FOR TH=S TO T
   LET X(TH)=XX+R*COS(TH)
   LET Y(TH)=YY-R*SIN(TH)
NEXT  TH
SET AREA COLOR C
MAT PLOT AREA : X,Y
END SUB

EXTERNAL FUNCTION CVI(A$)
OPTION ARITHMETIC NATIVE
OPTION CHARACTER BYTE
DECLARE NUMERIC A
LET A$=LEFT$(A$,2)
LET A=ORD(A$(1:1))+ORD(A$(2:2))*256
IF A>32767 THEN LET A=A-65536
LET CVI=A
END FUNCTION

EXTERNAL FUNCTION CVL(A$)
OPTION ARITHMETIC NATIVE
OPTION CHARACTER BYTE
DECLARE NUMERIC A
LET A$=LEFT$(A$,4)
LET A=ORD(A$(1:1))+ORD(A$(2:2))*256+ORD(A$(3:3))*256^2+ORD(A$(4:4))*256^3
IF A>=2^31-1 THEN LET A=A-2^32
LET CVL=A
END FUNCTION

EXTERNAL SUB CDFT(N,WR,WI,AR(),AI())
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM A(N)
FOR I=0 TO N/2-1
   LET A(2*I)=AR(I)
   LET A(2*I+1)=AI(I)
NEXT I
LET WMR=WR
LET WMI=WI
LET M=N
DO WHILE M>4
   LET L=M/2
   LET WKR=1
   LET WKI=0
   LET WDR=1-2*WMI*WMI
   LET WDI=2*WMI*WMR
   LET SS=2*WDI
   LET WMR=WDR
   LET WMI=WDI
   FOR J=0 TO N-M STEP M
      LET I=J+L
      LET XR=A(J)-A(I)
      LET XI=A(J+1)-A(I+1)
      LET A(J)=A(J)+A(I)
      LET A(J+1)=A(J+1)+A(I+1)
      LET A(I)=XR
      LET A(I+1)=XI
      LET XR=A(J+2)-A(I+2)
      LET XI=A(J+3)-A(I+3)
      LET A(J+2)=A(J+2)+A(I+2)
      LET A(J+3)=A(J+3)+A(I+3)
      LET A(I+2)=WDR*XR-WDI*XI
      LET A(I+3)=WDR*XI+WDI*XR
   NEXT J
   FOR K=4 TO L-4 STEP 4
      LET WKR=WKR-SS*WDI
      LET WKI=WKI+SS*WDR
      LET WDR=WDR-SS*WKI
      LET WDI=WDI+SS*WKR
      FOR J=K TO N-M+K STEP M
         LET I=J+L
         LET XR=A(J)-A(I)
         LET XI=A(J+1)-A(I+1)
         LET A(J)=A(J)+A(I)
         LET A(J+1)=A(J+1)+A(I+1)
         LET A(I)=WKR*XR-WKI*XI
         LET A(I+1)=WKR*XI+WKI*XR
         LET XR=A(J+2)-A(I+2)
         LET XI=A(J+3)-A(I+3)
         LET A(J+2)=A(J+2)+A(I+2)
         LET A(J+3)=A(J+3)+A(I+3)
         LET A(I+2)=WDR*XR-WDI*XI
         LET A(I+3)=WDR*XI+WDI*XR
      NEXT J
   NEXT K
   LET M=L
LOOP
IF M>2 THEN
   FOR J=0 TO N-4 STEP 4
      LET XR=A(J)-A(J+2)
      LET XI=A(J+1)-A(J+3)
      LET A(J)=A(J)+A(J+2)
      LET A(J+1)=A(J+1)+A(J+3)
      LET A(J+2)=XR
      LET A(J+3)=XI
   NEXT J
END IF
IF N>4  THEN CALL BITRV2(N,A)
FOR I=0 TO N/2-1
   LET AR(I)=A(2*I)/SQR(N)
   LET AI(I)=A(2*I+1)/SQR(N)
NEXT I
END SUB

EXTERNAL SUB BITRV2(N,A())
OPTION ARITHMETIC NATIVE
LET M=N/4
LET M2=2*M
LET N2=N-2
LET K=0
FOR J=0 TO M2-4 STEP 4
   IF J<K THEN
      LET XR=A(J)
      LET XI=A(J+1)
      LET A(J)=A(K)
      LET A(J+1)=A(K+1)
      LET A(K)=XR
      LET A(K+1)=XI
   ELSEIF J>K THEN
      LET J1=N2-J
      LET K1=N2-K
      LET XR=A(J1)
      LET XI=A(J1+1)
      LET A(J1)=A(K1)
      LET A(J1+1)=A(K1+1)
      LET A(K1)=XR
      LET A(K1+1)=XI
   END IF
   LET K1=M2+K
   LET XR=A(J+2)
   LET XI=A(J+3)
   LET A(J+2)=A(K1)
   LET A(J+3)=A(K1+1)
   LET A(K1)=XR
   LET A(K1+1)=XI
   LET L=M
   DO WHILE K>=L
      LET K=K-L
      LET L=L/2
   LOOP
   LET K=K+L
NEXT J
END SUB

EXTERNAL FUNCTION HANNING(X,N)
OPTION ARITHMETIC NATIVE
LET  HANNING=.5-.5*COS(2*PI*X/N) !'ハニング窓関数
END FUNCTION
 

オーディオスペクトラム画像

 投稿者:しばっち  投稿日:2018年 5月20日(日)20時13分26秒
返信・引用
  You Tubeを見ているとオーディオスペクトラム動画なるものを見かけることがあります。
これらの動画は専用のソフトや市販ソフトを使って作られているようです。
専用のソフト等には到底及びませんが、十進BASICでオーディオスペクトラム画像を書出します。
ここではGIFアニメではなく動画ファイル(mp4形式)として作成します。(画像サンプルを参照)

その補助ツール(メインツール!?)にffmpegを使用します。
ffmpeg.exeはコマンドラインアプリです。

https://www.ffmpeg.org

FFmpeg Buildsからwindows版 version 4.0 windows 64bit staticをダウンロードしました。


まず、WAVファイルを用意します。
WAVファイルをそのままで保存している方は少ないと思いますので
ないなら変換してWAVファイルを作成します。
想定しているWAVファイルは44.1kHzサンプリング 16bit 2チャンネルです。

この変換にffmpegが使用できます。コマンドプロンプトを起動して

mp3やflac,oggなどの音楽ファイルから

ffmpeg -i sample.mp3 sample.wav
ffmpeg -i sample.mp3 -acodec pcm_s16le -ar 44100 -ac 2 sample.wav (16bit pcm, 44.1kHz, 2channel)

mp4やmpg,flvなどの動画ファイルから

ffmpeg -i sample.flv -vn sample.wav (-vn videoなし)

ffmpeg --help とすると、簡単なヘルプが表示されます。


曲長として 4~5分程度を想定しています。ドライブ残量に気を付けてください。
5分の曲とすると

44100Hz * 16bit/8 * 2ch * 300秒 = 52920000 byte
フレームレートは30fpsを想定しています。300秒 * 30fps =9000枚の画像(jpg)を十進BASICで書き出します。

用意ができたら十進BASICを実行して、WAVファイルを読み込みます。
但し、全てのWAVファイルに対応しているわけではありません。
実行には曲長の数倍の時間がかかります。


出来上がった5桁の連番画像ファイルと音楽ファイルで、mp4ファイルに変換します。
フレームレート30fps ビデオコーデックh264 オーディオコーデックaac を指定します。

ffmpeg -r 30 -i image%05d.jpg -i sample.wav -vcodec h264 -acodec aac sample.mp4
ffmpeg -framerate 30 -i image%05d.jpg -i sample.wav -r 30 -vcodec libx264 -acodec aac sample.mp4

また、設定しだいでmpg,flvやwebmなどが作成できます。

大量に画像を作成しますので動作をテストモードにしています。
変数TESTを書き換えてください。


実際に動画作成してみました。
曲長 237.0989569161秒 画像 7112枚を出力 画像1枚あたり多くて100kB前後 少なくて20kB前後
mp4に変換すると、およそ30MB程のmp4動画ができました。
                                              (チャン! チャン!)

OPTION ARITHMETIC NATIVE
OPTION BASE 0
OPTION CHARACTER BYTE
DIM A$(20)
LET XSIZE=800    !'画像サイズ
LET YSIZE=400
LET FRAMERATE=30 !'動画フレームレート 30fps
LET DISPLAYMODE=0
!'INPUT  PROMPT "DISPLAY MODE (0 - 9)=":DISPLAYMODE
LET TEST=0        !'画像の書出し 0以外にすると画像書出し。
LET SCALE=YSIZE/XSIZE
SET TEXT BACKGROUND "OPAQUE"
FILE GETNAME F$,"WAVファイル|*.WAV"
IF F$="" THEN STOP
FILE SPLITNAME (F$) PATH$, NAME$, EXT$
OPEN #1:NAME F$,ACCESS INPUT !'wavファイル読み込み
FOR I=1 TO 12
   CHARACTER INPUT #1:A$(I)
NEXT I
IF A$(1)&A$(2)&A$(3)&A$(4)<>"RIFF" THEN
   PRINT "WAVファイルではありません"
   CLOSE #1
   STOP
END IF
LET WAVEFILESIZE=CVL(A$(5)&A$(6)&A$(7)&A$(8))
IF A$(9)&A$(10)&A$(11)&A$(12)<>"WAVE" THEN
   PRINT "WAVファイルではありません"
   CLOSE #1
   STOP
END IF
DO
   FOR I=1 TO 4
      CHARACTER INPUT #1:A$(I)
   NEXT I
   SELECT CASE A$(1)&A$(2)&A$(3)&A$(4)
   CASE "fmt "
      FOR I=1 TO 4
         CHARACTER INPUT #1:A$(I)
      NEXT I
      LET HEADERSIZE=CVL(A$(1)&A$(2)&A$(3)&A$(4))
      FOR I=1 TO HEADERSIZE
         CHARACTER INPUT #1:A$(I)
      NEXT I
      LET WAVETYPE=CVI(A$(1)&A$(2))
      IF WAVETYPE<>1 THEN
         PRINT "対応していません"
         CLOSE #1
         STOP
      END IF
      LET CHANNEL=CVI(A$(3)&A$(4))
      LET SAMPLINGFREQ=CVL(A$(5)&A$(6)&A$(7)&A$(8))
      LET DATARATE=CVL(A$(9)&A$(10)&A$(11)&A$(12))
      LET SAMPLESIZE=CVI(A$(13)&A$(14))
      LET SAMPLEBIT=CVI(A$(15)&A$(16))
   CASE "data"
      FOR I=1 TO 4
         CHARACTER INPUT #1:A$(I)
      NEXT I
      LET PCMSIZE=CVL(A$(1)&A$(2)&A$(3)&A$(4))
      LET SECOND=PCMSIZE/DATARATE
      !' LET NUM=INT(SAMPLINGFREQ*SECOND)
      LET NUM=PCMSIZE/SAMPLESIZE
      EXIT DO
   CASE "fact"
      FOR I=1 TO 8
         CHARACTER INPUT #1:A$(I)
      NEXT I
      !' LET SIZE=CVL(A$(1)&A$(2)&A$(3)&A$(4))
      !' LET SAMPLE=CVL(A$(5)&A$(6)&A$(7)&A$(8))
   CASE "LIST"
      FOR I=1 TO 4
         CHARACTER INPUT #1:A$(I)
      NEXT I
      LET SIZE=CVL(A$(1)&A$(2)&A$(3)&A$(4))
      FOR I=1 TO SIZE
         CHARACTER INPUT #1:DMY$ !'空読み
      NEXT I
   CASE ELSE
      PRINT "対応していません"
      CLOSE #1
      STOP
   END SELECT
LOOP
PRINT "ファイルサイズ ";WAVEFILESIZE;"Byte" !' データ表示
PRINT "WAVEタイプ ";WAVETYPE
PRINT "チャンネル数 ";CHANNEL;"ch"
PRINT "サンプリング周波数 ";SAMPLINGFREQ;"Hz"
PRINT "データレート ";DATARATE;"Byte"
PRINT "ブロックサイズ ";SAMPLESIZE;"Byte"
PRINT "サンプリングビット ";SAMPLEBIT;"bit"
PRINT "PCMサイズ ";PCMSIZE;"Byte"
PRINT "データ数 ";PCMSIZE/SAMPLESIZE;INT(SAMPLINGFREQ*SECOND)
PRINT "時間 ";SECOND;"秒"
PRINT "描画モード ";DISPLAYMODE
PRINT "出力画像 ";INT(SECOND*FRAMERATE);"枚"
SET POINT STYLE 1
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0
SET COLOR MIX(1) 0,0,1
SET COLOR MIX(2) 1,0,0
SET COLOR MIX(3) 1,0,1
SET COLOR MIX(4) 0,1,0
SET COLOR MIX(5) 0,1,1
SET COLOR MIX(6) 1,1,0
SET COLOR MIX(7) 1,1,1
SET BITMAP SIZE XSIZE,YSIZE
LET N=INT(SAMPLINGFREQ/FRAMERATE) !'1フレームあたりのポイント数
LET NN=2^INT(LOG2(N)+1) !'FFTポイント数
DIM RMAX(NN),LR(NN),RR(NN),FR(NN),FI(NN)
LET B$="  "
FOR K=0 TO NUM-1
   FOR CH=1 TO CHANNEL
      FOR J=1 TO SAMPLEBIT/8
         CHARACTER INPUT #1:B$(J:J)
      NEXT J
      IF SAMPLEBIT=8 THEN LET DAT=ORD(B$(1:1))-128
      IF SAMPLEBIT=16 THEN LET DAT=CVI(B$)
      LET J=MOD(K,N)
      SELECT CASE CH
      CASE 1
         LET LR(J)=DAT
      CASE 2
         LET RR(J)=DAT
      CASE ELSE
      END SELECT
   NEXT CH
   IF CHANNEL=1 THEN MAT RR=LR !'モノラルならコピー
   IF K>0 AND J=0 THEN !'1フレーム描画 30フレームで1秒分
      LET COUNT=COUNT+1
      SET DRAW MODE HIDDEN !'ちらつき防止
      SELECT CASE DISPLAYMODE !'画面表示
      CASE 0
         CALL DISPLAY0
      CASE 1
         CALL DISPLAY1
      CASE 2
         CALL DISPLAY2
      CASE 3
         CALL DISPLAY3
      CASE 4
         CALL DISPLAY4
      CASE 5
         CALL DISPLAY5
      CASE 6
         CALL DISPLAY6
      CASE 7
         CALL DISPLAY7
      CASE 8
         CALL DISPLAY8
      CASE 9
         CALL DISPLAY9
      END SELECT
      SET DRAW MODE EXPLICIT
      IF TEST<>0 THEN GSAVE PATH$&"image"&RIGHT$("00000"&STR$(COUNT),5)&".jpg" !'画像書出し
   END IF
NEXT K
CLOSE #1
PRINT "ffmpeg -r";FRAMERATE;"-i ";CHR$(34);PATH$;"image%05d.jpg";CHR$(34);" -i ";CHR$(34);F$;CHR$(34);" -vcodec h264 -b:v 1000k -acodec aac -b:a 128k ";CHR$(34);PATH$;NAME$;".mp4";CHR$(34)

SUB DISPLAY1
   CALL GCLEAR
   CALL TIMER(.9,1,0,.05,2,12) !'タイマー
   CALL WAVE(RR,.51,1,.05,.95,5) !'右チャンネル(右)
   CALL WAVE(LR,0,.49,.05,.95,5) !'左チャンネル(左)
END SUB

SUB DISPLAY2
   CALL GCLEAR
   CALL TIMER(.9,1,0,.05,2,12) !'タイマー
   CALL WAVE2(LR,0,1,.6,.9,4) !'左チャンネル(上)
   CALL WAVE2(RR,0,1,.1,.4,4) !'右チャンネル(下)
END SUB

SUB DISPLAY3
   CALL GCLEAR
   CALL TIMER(.9,1,0,.05,2,12)
   CALL FREQ(LR,0,1,.6,.9,5,2) !'左チャンネル(上)
   CALL FREQ(RR,0,1,.1,.4,5,2) !'右チャンネル(下)
END SUB

SUB DISPLAY4
   CALL GCLEAR
   CALL POWER(LR,.05,.48,.05,.95,4,2,8,"MODE2","RIGHT") !'左チャンネル(左)
   CALL POWER(RR,.52,.95,.05,.95,4,2,8,"MODE2","LEFT") !'右チャンネル(右)
END SUB

SUB DISPLAY5
   CALL GCLEAR
   CALL WAVE3(LR,0,1,.6,.9,6,2) !'左チャンネル(上)
   CALL WAVE3(RR,0,1,.1,.4,6,2) !'右チャンネル(下)
END SUB

SUB DISPLAY6
   CALL GCLEAR
   CALL VOLUMEBAR(LR,.1,.9,.7,.8,4,2,"MODE2","LEFT") !'左チャンネル(上)
   CALL VOLUMEBAR(RR,.1,.9,.2,.3,4,2,"MODE2","LEFT") !'右チャンネル(下)
END SUB

SUB DISPLAY7
   CALL GCLEAR
   CALL VOLUMECIRCLE(LR,.05,.45,.1,.9,4,2) !'左チャンネル(左)
   CALL VOLUMECIRCLE(RR,.55,.95,.1,.9,4,2) !'右チャンネル(右)
END SUB

SUB DISPLAY8
   CALL GCLEAR
   CALL XYPLOT(RR,.51,1,.05,.95,7) !'右チャンネル(右)
   CALL XYPLOT(LR,0,.49,.05,.95,7) !'左チャンネル(左)
END SUB

SUB DISPLAY9
   CALL GCLEAR
   CALL POWERDISK(RR,.51,1,.05,.95,5,2,8,"MODE2") !'右チャンネル(右)
   CALL POWERDISK(LR,0,.49,.05,.95,5,2,8,"MODE2") !'左チャンネル(左)
END SUB

SUB DISPLAY0
   CALL GCLEAR
   CALL TIMER(.9,1,0,.04,2,10)
   CALL TIMEBAR(.1,.8,0,.04,3)
   CALL VOLUMEBAR(LR,.02,.49,.05,.08,1,2,"MODE2","LEFT")
   CALL POWER(LR,.02,.49,.1,.5,4,2,8,"MODE2","BOTTOM")
   CALL FREQ(LR,.02,.49,.52,.75,6,2)
   CALL WAVE(LR,.02,.49,.76,.99,5)
   CALL VOLUMEBAR(RR,.51,.98,.05,.08,1,2,"MODE2","LEFT")
   CALL POWER(RR,.51,.98,.1,.5,4,2,8,"MODE2","BOTTOM")
   CALL FREQ(RR,.51,.98,.52,.75,6,2)
   CALL WAVE(RR,.51,.98,.76,.99,5)
END SUB

SUB GCLEAR
   SET VIEWPORT 0,1,0,SCALE
   SET WINDOW 0,1,0,SCALE
   SET AREA COLOR 7
   PLOT AREA:0,0;1,0;1,SCALE;0,SCALE
END SUB

SUB TIMER(XS,XE,YS,YE,COL,HEIGHT)
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE !'VIEWPORTの領域を長方形でも 左下(0,0) 右上(1,1)とする
   SET WINDOW XS,XE,YS,YE
   SET TEXT FONT "",HEIGHT
   SET TEXT COLOR COL
   PLOT TEXT ,AT XS,YS:RIGHT$("0"&STR$(INT(COUNT/FRAMERATE/60)),2)&":"&RIGHT$("0"&STR$(MOD(INT(COUNT/FRAMERATE),60)),2)&"."&RIGHT$("0"&STR$(MOD(COUNT,FRAMERATE)),2)
END SUB

SUB TIMEBAR(XS,XE,YS,YE,COL)
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET XS=0
   LET XE=NUM-1
   LET YS=0
   LET YE=1
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   CALL BOXFULL(0,0,K,1,COL)
END SUB

SUB WAVE(RR(),XS,XE,YS,YE,COL) !'波形
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET XS=0
   LET XE=N-1
   LET YS=-2^(SAMPLEBIT-1)
   LET YE=-YS
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   SET LINE COLOR COL
   FOR I=0 TO N-1
      PLOT LINES:I,RR(I);
   NEXT I
END SUB

SUB WAVE2(RR(),XS,XE,YS,YE,COL)
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET XS=0
   LET XE=N-1
   LET YS=-2^(SAMPLEBIT-1)
   LET YE=-YS
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   SET LINE COLOR COL
   FOR I=0 TO N-1
      PLOT LINES:I,-RR(I);I,RR(I)
   NEXT I
END SUB

SUB WAVE3(RR(),XS,XE,YS,YE,COL,MAXCOL)
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET XS=0
   LET XE=N-1
   LET YE=2^(SAMPLEBIT-1)
   LET YS=0
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   IF COUNT=1 OR MOD(COUNT,INT(FRAMERATE/2))=0 THEN !'フレームレートの半分で最大値クリア
      MAT RMAX=ZER
   END IF
   SET LINE COLOR MAXCOL
   FOR I=0 TO N-1
      PLOT LINES:I,0;I,RMAX(I)
   NEXT I
   PLOT LINES
   SET LINE COLOR COL
   FOR I=0 TO N-1
      LET R=ABS(RR(I))
      LET RMAX(I)=MAX(RMAX(I),R)
      PLOT LINES:I,0;I,R
   NEXT I
END SUB

SUB FREQ(RR(),XS,XE,YS,YE,COL,MAXCOL) !'周波数
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET XS=0
   LET XE=NN/2-1
   LET YS=0
   LET YE=20*LOG10(2^(SAMPLEBIT-1))
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   IF COUNT=1 OR MOD(COUNT,INT(FRAMERATE/2))=0 THEN
      MAT RMAX=ZER
   END IF
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN) !'ハニング窓関数を掛ける
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI) !'FFT
   SET LINE COLOR COL
   FOR I=0 TO NN/2-1
      LET R=10*LOG10(FR(I)^2+FI(I)^2+1)
      LET RMAX(I)=MAX(RMAX(I),R)
      PLOT LINES:I,R;
   NEXT I
   PLOT LINES
   SET LINE COLOR MAXCOL
   FOR I=0 TO NN/2-1
      PLOT LINES:I,RMAX(I); !'最大値描画
   NEXT I
END SUB

SUB POWER(RR(),XS,XE,YS,YE,COL,MAXCOL,LEV,MODE$,DIRECTION$) !'スペクトル
   SET VIEWPORT XS,XE,YS*SCALE,YE*SCALE
   LET P=10*LOG10(2^(SAMPLEBIT-1))
   SELECT CASE DIRECTION$
   CASE "BOTTOM"
      LET YS=0
      LET YE=P
      LET XE=NN/2-1
      LET XS=0
   CASE "TOP"
      LET YS=P
      LET YE=0
      LET XS=0
      LET XE=NN/2-1
   CASE "LEFT"
      LET XS=0
      LET XE=P
      LET YS=0
      LET YE=NN/2-1
   CASE "RIGHT"
      LET XS=P
      LET XE=0
      LET YS=0
      LET YE=NN/2-1
   END SELECT
   SET WINDOW XS,XE,YS,YE
   CALL BOXFULL(XS,YS,XE,YE,0)
   FOR I=0 TO NN-1
      LET FR(I)=RR(I)*HANNING(I,NN)
   NEXT I
   MAT FI=ZER
   CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
   LET T=NN/LEV/2
   FOR I=0 TO NN/2-1 STEP T
      LET R=0
      FOR P=0 TO T-1
         LET R=R+10*LOG10(FR(I+P)^2+FI(I+P)^2+1)
      NEXT P
      LET R=R/T !'平均値
      SELECT CASE DIRECTION$
      CASE "TOP","BOTTOM"
         SELECT CASE MODE$
         CASE "MODE1"
            IF R>MAX(YS,YE)/32 THEN
               FOR Y=0 TO R STEP MAX(YS,YE)/32
                  IF Y<MAX(YS,YE)*.8 THEN CALL BOXFULL(I,Y,T*3/4+I,Y+MAX(YS,YE)/32,COL) ELSE CALL BOXFULL(I,Y,T*3/4+I,Y+MAX(YS,YE)/32,MAXCOL)
               NEXT Y
            END IF
         CASE "MODE2"
            IF R>MAX(YS,YE)/16 THEN
               FOR Y=0 TO R STEP MAX(YS,YE)/16
                  IF Y<MAX(YS,YE)*.8 THEN CALL BOXFULL(I,Y,T*3/4+I,Y+MAX(YS,YE)/16*.8,COL) ELSE CALL BOXFULL(I,Y,T*3/4+I,Y+MAX(YS,YE)/16*.8,MAXCOL)
               NEXT Y
            END IF
         END SELECT
      CASE "LEFT","RIGHT"
         SELECT CASE MODE$
         CASE "MODE1"
            IF R>MAX(XS,XE)/32 THEN
               FOR X=0 TO R STEP MAX(XS,XE)/32
                  IF X<MAX(XS,XE)*.8 THEN CALL BOXFULL(X,I,X+MAX(XS,XE)/32,T*3/4+I,COL) ELSE CALL BOXFULL(X,I,X+MAX(XS,XE)/32,T*3/4+I,MAXCOL)
               NEXT X
            END IF
         CASE "MODE2"
            IF R>MAX(XS,XE)/16 THEN
               FOR X=0 TO R STEP MAX(XS,XE)/16
                  IF X<MAX(XS,XE)*.8 THEN CALL BOXFULL(X,I,X+MAX(XS,XE)/16*.8,T*3/4+I,COL) ELSE CALL BOXFULL(X,I,X+MAX(XS,XE)/16*.8,T*3/4+I,MAXCOL)
               NEXT  X
            END IF
         END SELECT
      END SELECT
   NEXT I
END SUB


続く
 

[025] イオン結晶の高速化-分子動力学法プログラム

 投稿者:mike  投稿日:2018年 5月12日(土)11時02分32秒
返信・引用
  イオン結晶の高速化-分子動力学法プログラム025fasterIonMD2D.basを公開します。([023]の高速化バージョンです。)
本プログラムは十進BASIC 8.0.0 / macOS 10.7.5 でテストしました。

! ========= molecular dynamics 2D - ion ==========
!
! 025NaClIMD2D.bas
!   Copyright(C) 2018 Mitsuru Ikeuchi
!   Released under the MIT license ( https://opensource.org/licenses/MIT )
!
! ver 0.0.1   2018.05.11  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB imd2d.setInitialCondition, imd2d.moveParticles, imd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET ddTemp = 0      !contTemp = contTemp + ddTemp
LET drawMode = 0    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
LET pauseFlag = 0   !if pauseFlag=0, CALL moveParticles(tempMode,contTemp)

!setInitialCondition(material,boxSizeInNM,contTemp)
CALL setInitialCondition(1,12.0,contTemp) !material = 1:NaCl 2:MgO 3:CaO 4:BaO 5:NaF 6:KF 7:KCl

DO
   LET contTemp = contTemp + ddTemp
   IF contTemp>3000 THEN LET contTemp = 3000
   IF contTemp<10 THEN LET contTemP = 10
   IF pauseFlag=0 THEN CALL moveParticles(tempMode,contTemp) ELSE WAIT DELAY 0.05
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" OR S$="7" THEN
      LET material = VAL(S$)
      LET tempMode = 0
      LET contTemp = 300
      CALL setInitialCondition(material,12.0,contTemp)
   ELSEIF S$="T" OR S$="t" THEN
      LET tempMode = MOD(tempMode+1,2)
      IF tempMode=0 THEN LET ddTemp = 0
   ELSEIF S$="K" OR S$="k" THEN
      LET tempMode = 1
      IF ddtemp=0 THEN LET ddTemp = 1 ELSE LET ddTemp = 0
   ELSEIF S$="J" OR S$="j" THEN
      LET tempMode = 1
      IF ddTemp=0 THEN LET ddTemp = -1 ELSE LET ddTemp = 0
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$=" " THEN
      LET pauseFlag = MOD(pauseFlag+1,2)
   END IF
LOOP UNTIL S$=CHR$(27)
END

EXTERNAL FUNCTION INKEY$ !--- from decimal BASIC library inkey$.bas
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- molecular dynamics 2D - ion ----------
!
!  method: velocity Verlet ( F=m*d^2r/dt^2 -> r(t+dt)=r(t)+v*dt,v=v+(F/m)*dt )
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!
!  force: ion f(r) = fc + fr + fa
!        fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force
!        fr = 6.9742e-11*EXP((a-r)/b) !repulsive force
!        fa = -6.9742e-21*(c/r^6) !attractive force
!
MODULE imd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC ionKind1,ionKind2, sysTime, dt, nMolec, xMax, yMax, Nsx,Nsy, rCutoff, hh
SHARE NUMERIC xx(5000),yy(5000)             ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(5000),vy(5000)             ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(5000),ffy(5000)           ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(5000),mas(5000)          ! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(5000,0 TO 100)            ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC section(0 TO 101,0 TO 101,0 TO 20) !use pre-registration
SHARE NUMERIC ion_m(0 TO 17),ion_z(0 TO 17) ! ion mass,ion charge
SHARE NUMERIC ion_a(0 TO 17),ion_b(0 TO 17) ! ion force param a,ion force param b
SHARE NUMERIC ion_c(0 TO 17),ion_r(0 TO 17) ! ion force param c, ion radius
SHARE NUMERIC ion_col(0 TO 17)              ! ion draw color
SHARE STRING ion_str$(0 TO 17)              ! ion string
SHARE NUMERIC forceTable(0 TO 17, 0 TO 17,0 TO 1001) ! force table

LET ionKind1 = 3           ! ion kind1 - 3:Na+
LET ionKind2 = 7           ! ion kind2 - 7:Cl-
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 2.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET Nsx = 100
LET Nsy = 100
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET ion_z(0)=0.0           ! if ion_z(0)=0.0 then set ion data

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(material,boxSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setIonData,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setNaClTypeBlock
   RANDOMIZE
   LET sysTime = 0.0
   CALL setIonData
   CALL setForceTable
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   IF material=1 THEN !NaCl
      LET ionKind1 = 3 !Na+
      LET ionKind2 = 7 !Cl-
      LET lattice = 5.6407e-10
   ELSEIF material=2 THEN !MgO
      LET ionKind1 = 4 !Mg++
      LET ionKind2 = 1 !O--
      LET lattice = 4.212e-10*0.94 !0.94: correction factor
   ELSEIF material=3 THEN !CaO
      LET ionKind1 = 9 !Ca++
      LET ionKind2 = 1 !O--
      LET lattice = 4.80e-10*0.94 !0.94: correction factor
   ELSEIF material=4 THEN !BaO
      LET ionKind1 = 12 !Ba++
      LET ionKind2 = 1 !O--
      LET lattice = 5.536e-10
   ELSEIF material=5 THEN !NaF
      LET ionKind1 = 3 !Na+
      LET ionKind2 = 2 !F-
      LET lattice = 4.62e-10
   ELSEIF material=6 THEN !KF
      LET ionKind1 = 8 !K+
      LET ionKind2 = 2 !F-
      LET lattice = 5.34e-10
   ELSEIF material=7 THEN !KCl
      LET ionKind1 = 8 !K+
      LET ionKind2 = 7 !Cl-
      LET lattice = 6.29e-10
   END IF
   LET nL = int(0.8*xMax/lattice)
   LET s = 0.5*(xMax - lattice*nL)
   !            setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos)
   LET nMolec = setNaClTypeBlock(1, ionKind1, ionKind2, nL, nL, lattice, s, s)
   CALL ajustVelocity(contTemp)
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setIonData
! ion potential data
!      0 mass,1 charge,  2  a ,     3  b , 4 c , 5  r-ion, 6 color 7 str$
   DATA 10.81,  3.0, 0.720e-10, 0.080e-10,  0.0, 0.23e-10, 13 , "B+++"   !0
   DATA 16.00, -2.0, 1.626e-10, 0.085e-10, 20.0, 1.40e-10,  2 , "O--"    !1
   DATA 19.00, -1.0, 1.565e-10, 0.085e-10, 20.0, 1.33e-10,  2 , "F-"     !2
   DATA 22.99,  1.0, 1.260e-10, 0.080e-10, 20.0, 1.02e-10,  4 , "Na+"    !3
   DATA 24.31,  2.0, 1.161e-10, 0.080e-10, 10.0, 0.72e-10,  4 , "Mg++"   !4
   DATA 26.98,  3.0, 1.064e-10, 0.080e-10,  2.0, 0.53e-10, 13 , "Al+++"  !5
   DATA 28.09,  4.0, 1.012e-10, 0.080e-10,  0.0, 0.40e-10,  7 , "Si++++" !6
   DATA 35.45, -1.0, 1.950e-10, 0.090e-10, 30.0, 1.81e-10,  2 , "Cl-"    !7
   DATA 39.10,  1.0, 1.595e-10, 0.080e-10, 15.0, 1.38e-10,  4 , "K+"     !8
   DATA 40.08,  2.0, 1.414e-10, 0.080e-10, 10.0, 1.00e-10,  4 , "Ca++"   !9
   DATA 47.88,  4.0, 1.235e-10, 0.080e-10,  0.0, 0.61e-10,  7 , "Ti++++" !10
   DATA 87.62,  2.0, 1.632e-10, 0.080e-10, 15.0, 1.16e-10,  4 , "Sr++"   !11
   DATA 137.3,  2.0, 1.820e-10, 0.080e-10, 20.0, 1.36e-10,  4 , "Ba++"   !12
   DATA 4.003,  0.0, 1.200e-10, 0.110e-10, 4.76, 1.28e-10,  3 , "He"     !13
   DATA 20.18,  0.0, 1.415e-10, 0.112e-10,11.03, 1.37e-10,  3 , "Ne"     !14
   DATA 39.95,  0.0, 1.878e-10, 0.117e-10,38.53, 1.70e-10,  3 , "Ar"     !15
   DATA 83.80,  0.0, 2.041e-10, 0.130e-10,55.33, 1.83e-10,  3 , "Kr"     !16
   DATA 131.3,  0.0, 2.258e-10, 0.145e-10,85.55, 1.99e-10,  3 , "Xe"     !17

   IF ion_z(0)=0.0 THEN
      FOR i=0 TO 17
         READ ion_m(i),ion_z(i),ion_a(i),ion_b(i),ion_c(i),ion_r(i),ion_col(i),ion_str$(i)
         LET ion_m(i) = ion_m(i)*1.661e-27
         LET ion_z(i) = ion_z(i)*1.602e-19
      NEXT i
   END IF
END SUB

EXTERNAL FUNCTION setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos)
   DECLARE EXTERNAL SUB setParticle
   LET ipp = ii
   LET a = lattice/2.0
   FOR i=0 TO 2*nx-1
      FOR j=0 TO 2*ny-1
         LET x = xPos + a*i
         LET y = yPos + a*j
         IF MOD(i+j,2)=0 THEN
            LET knd = knd1
         ELSE
            LET knd = knd2
         END IF
         CALL setParticle(ipp, knd, x, y)
         LET ipp = ipp + 1
      NEXT j
   NEXT i
   LET setNaClTypeBlock = ipp - 1
END FUNCTION

EXTERNAL SUB setParticle(i, knd, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
   LET mas(i) = ion_m(knd)
   LET kind(i) = knd
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   LET eForceConst = 1.0/(4.0*PI*8.8542e-12) !epsilon0=8.8542e-12
   FOR ki=0 TO 17
      FOR kj=0 TO 17
         LET a = ion_a(ki)+ion_a(kj)
         LET b = ion_b(ki)+ion_b(kj)
         LET c = ion_c(ki)*ion_c(kj)*1.0e-60
         LET zi = ion_z(ki)
         LET zj = ion_z(kj)
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force
            LET fr = 6.9742e-11*EXP((a-r)/b) !repulsive force
            LET fa = -6.9742e-21*(c/(r*r*r*r*r*r)) !attractive force
            LET forceTable(ki,kj,ir) = cutoff(r)*(fc + fr + fa)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec,registration
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registration
   !CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

EXTERNAL SUB registration !faster registration
   DECLARE EXTERNAL SUB preRegistration
   CALL preRegistration
   LET rreg = rCutoff+20*2000*dt
   LET rreg2 = rreg*rreg
   FOR ipp=1 TO nMolec-1
      LET kp = 1
      LET i0 = INT(Nsx*(xx(ipp)-rreg)/xMax)
      IF (i0<0) THEN LET i0 = 0
      LET i1 = INT(Nsx*(xx(ipp)+rreg)/xMax )
      IF (i1>=Nsx) THEN LET i1 = Nsx-1
      LET j0 = INT(Nsy*(yy(ipp)-rreg)/yMax )
      IF (j0<0) THEN LET j0 = 0
      LET j1 = INT(Nsy*(yy(ipp)+rreg)/yMax )
      IF (j1>=Nsy) THEN LET j1 = Nsy-1
      FOR i=i0 TO i1
         FOR j=j0 TO j1
            FOR iq=1 TO section(i,j,0)
               LET jp = section(i,j,iq)
               IF (jp>ipp) THEN
                  LET r2=(xx(ipp)-xx(jp))*(xx(ipp)-xx(jp))+(yy(ipp)-yy(jp))*(yy(ipp)-yy(jp))
                  IF (r2<rreg2) THEN
                     LET reg(ipp,kp) = jp
                     LET kp = kp + 1
                  END IF
               END IF
            NEXT iq
         NEXT j
      NEXT i
      LET reg(ipp,0) = kp
   NEXT ipp
END SUB

EXTERNAL SUB preRegistration
   FOR i=0 TO Nsx-1
      FOR j=0 TO Nsy-1
         LET section(i,j,0) = 0
      NEXT j
   NEXT i
   FOR ipp=1 TO nMolec
      LET i = INT(Nsx*xx(ipp)/xMax)
      IF i>=Nsx THEN LET i = Nsx-1
      LET j = INT(Nsy*yy(ipp)/yMax)
      IF j>=Nsy THEN LET j = Nsy-1
      LET iq = section(i,j,0) + 1
      LET section(i,j,0) = iq
      LET section(i,j,iq) = ipp
   NEXT ipp
END SUB

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF

   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 :ion_str$(ionKind1)&","&ion_str$(ionKind2)
   PLOT TEXT, AT 100, 55 ,USING "  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional ion - molecular dynamics"
   PLOT TEXT, AT  50,480 :"[esc] exit  [space]pause/go  [D]change draw mode"
   PLOT TEXT, AT  50,460 :"[1]NaCl [2]MgO [3]CaO [4]BaO [5]NaF [6]KF [7]KCl"
   PLOT TEXT, AT  50,440 :"[T] toggle 0:adiabatic/1:constant temperature"
   PLOT TEXT, AT  50,420 :"temp control -> [J]coolDown/stop  [k]heatUp/stop"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      SET LINE COLOR ion_col(kind(i))
      DRAW circle WITH SCALE(ion_r(kind(i))*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            LET r0 = (ion_r(kind(i))+ion_r(kind(j)))
            IF r<1.1*r0 AND kind(i)<>kind(j) THEN
               SET LINE COLOR 8 !gray
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      IF kind(i)=ionKind1 THEN SET LINE COLOR 4 ELSE SET LINE COLOR 2  !4:red, 2:blue
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub

END MODULE

http://mike1336.web.fc2.com/

 

アップデート版「十進BASIC 第2掲示板のアーカイブス」

 投稿者:SECOND  投稿日:2018年 4月14日(土)18時37分37秒
返信・引用 編集済
  > No.4293[元記事へ]

サイズが、約47MB もあるので、ダウンロードの際、注意! 解凍すると、約97MB 。
http://neutro.la.coocan.jp/asm/Decimal-Basic-bbs2-180415.lzh

※Firefox で見ると 付属スレッドの左マージンが無くなっているのを、修正しました。
  Decimal-Basic-bbs2-180413.lzh → Decimal-Basic-bbs2-180415.lzh

■保存されている範囲。
 メインスレッド
  開始 No.1 新掲示板開設 投稿者:白石 和夫  投稿日:2008年 7月21日(月)09時38分46秒
      (
      )
  近況 No.4556 レイマーチング 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時14分28秒

 付属スレッド
  ◇スレッドが使えます(2)
  ◇Paract BASIC(21)
  ◇Amusement_Program(10)
  ◇改修予定のJIS非互換(3)
  ◇複数ページ長編プログラム(新規投稿)(16)
  ◇「十進BASIC第2掲示板」投稿記事リスト(17)
  ◇Full BASIC互換ライブラリ(8)
  ◇「十進BASIC掲示板過去ログ」インデックス(トピック)(17)
  ◇人の色覚の数理(14)
  ◇「十進BASIC掲示板過去ログ」インデックス(ツリー)(91)

 

レイマーチング

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時14分28秒
返信・引用
  ネット上のサンプルを基にレイマーチング法で球?をレンダリングしています。
二進モードで実行してください。

PUBLIC NUMERIC EPS,OFFSET
DIM POS(3),CPOS(3),CDIR(3),CUP(3),CSIDE(3),RAY(3),R1(3),R2(3),R3(3),LIGHTDIR(3)
DIM NORMAL(3),COL(3),C(3),V(3),R(3),CC(3),CAMERA(3)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
CALL VEC3(CDIR,0,0,1)
CALL VEC3(CUP,0,1,0)
CALL VEC3(CAMERA,0,1,0)
MAT CSIDE=CROSS(CDIR,CUP)
LET EPS=.01
LET TARGETDEPTH=1.3
LET OFFSET = EPS * 100
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      MAT R1=X*CSIDE
      MAT R2=Y*CUP
      MAT R3=TARGETDEPTH*CDIR
      MAT RAY=R1
      MAT RAY=RAY+R2
      MAT RAY=RAY+R3
      MAT CPOS=CAMERA
      CALL NORMALIZE(RAY)
      LET ALPHA=1
      MAT COL=ZER
      FOR I=0 TO 2
         CALL GETRAYCOLOR(CPOS, RAY, POS, NORMAL, HIT,CC)
         MAT C=ALPHA*CC
         MAT COL=COL+C
         LET ALPHA=ALPHA*.3
         CALL REFLECT(RAY, NORMAL ,R)
         CALL NORMALIZE(R)
         MAT RAY=R
         MAT CPOS=OFFSET*NORMAL
         MAT CPOS=CPOS+POS
         IF HIT=0 THEN EXIT FOR
      NEXT I
      CALL SETCOLOR(COL(1),COL(2),COL(3))
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB  ONREP(P(), INTERVAL,PP())
DIM Q(2)
LET Q(1)=MOD(P(1),INTERVAL)- INTERVAL * 0.5
LET Q(2)=MOD(P(3),INTERVAL)- INTERVAL * 0.5
CALL VEC3(PP,Q(1),P(2),Q(2))
END SUB

EXTERNAL  FUNCTION SPHEREDIST(P(),R)
DIM PP(3)
CALL ONREP(P,3,PP)
LET SPHEREDIST=LENGTH(PP)-R
END FUNCTION

EXTERNAL  FUNCTION FLOORDIST(P())
DIM V(3)
CALL VEC3(V,0,1,0)
LET FLOORDIST=DOT(P,V)+1
END FUNCTION

EXTERNAL  SUB MINVEC4(A(),B(),C())
IF A(4)<B(4) THEN MAT C=A ELSE MAT C=B
END SUB

EXTERNAL  FUNCTION CHECKEREDPATTERN(P())
LET U=1-IP(MOD(P(1),2))
LET V=1-IP(MOD(P(3),2))
IF U=1 AND V<1 OR U<1 AND V=1 THEN LET CHECKEREDPATTERN=.2 ELSE LET CHECKEREDPATTERN=1
END FUNCTION

EXTERNAL  SUB HSV2RGB(C(),RGB())
DIM K(4),P(3)
CALL VEC4(K,1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0)
FOR I=1 TO 3
   LET P(I)=ABS(FRACT(C(1)+K(I))*6-K(4))
NEXT I
FOR I=1 TO 3
   LET RGB(I)=C(3)*MIX(K(1),CLAMP(P(I)-K(1),0,1),C(2))
NEXT I
END SUB

EXTERNAL  FUNCTION SCENEDIST(P())
LET SCENEDIST=MIN(SPHEREDIST(P,1),FLOORDIST(P))
END FUNCTION

EXTERNAL  SUB SCENECOLOR(P(),PP())
DIM A(4),B(4),C(3),COL(3)
CALL VEC3(C,(P(3) + P(1)) / 9.0, 1.0, 1.0 )
CALL HSV2RGB(C,COL)
CALL VEC4(A,COL(1),COL(2),COL(3), SPHEREDIST(P,1.0))
LET L=CHECKEREDPATTERN(P)
CALL VEC4(B,.5*L,.5*L,.5*L,FLOORDIST(P))
CALL MINVEC4(A,B,PP)
END SUB

EXTERNAL  SUB GETNORMAL(P(),N())
DIM X1(3),X2(3),Y1(3),Y2(3),Z1(3),Z2(3),PX1(3),PX2(3),PY1(3),PY2(3),PZ1(3),PZ2(3)
CALL VEC3(X1,EPS,0,0)
CALL VEC3(X2,-EPS,0,0)
MAT PX1=P+X1
MAT PX2=P+X2
CALL VEC3(Y1,0,EPS,0)
CALL VEC3(Y2,0,-EPS,0)
MAT PY1=P+Y1
MAT PY2=P+Y2
CALL VEC3(Z1,0,0,EPS)
CALL VEC3(Z2,0,0,-EPS)
MAT PZ1=P+Z1
MAT PZ2=P+Z2
CALL VEC3(N,SCENEDIST(PX1)-SCENEDIST(PX2),SCENEDIST(PY1)-SCENEDIST(PY2),SCENEDIST(PZ1)-SCENEDIST(PZ2))
CALL NORMALIZE(N)
END SUB

EXTERNAL  FUNCTION GETSHADOW(RO(),RD())
DIM RAY(3),DMY(3)
LET R=1
LET SHADOWCOEF=.5
FOR T=0 TO 49
   MAT RAY=C*RD
   MAT RAY=RAY+RO
   LET H=SCENEDIST(RAY)
   IF H<EPS THEN
      LET GETSHADOW=SHADOWCOEF
      EXIT FUNCTION
   END IF
   IF C<>0 THEN LET R=MIN(R,H*16/C)
   LET C=C+H
NEXT T
LET GETSHADOW=1-SHADOWCOEF+R*SHADOWCOEF
END FUNCTION

EXTERNAL  SUB GETRAYCOLOR(ORIGIN(),RAY(),POS(),NORMAL(),HIT,COL())
DIM P(3),S(3),V(3),LIGHTDIR(3),CL(4)
MAT POS=ORIGIN
CALL VEC3(LIGHTDIR,-0.48666426339228763, 0.8111071056538127, -0.3244428422615251)
FOR I=0 TO 63
   LET DIST = SCENEDIST(POS)
   LET DEPTH =DEPTH+ DIST
   MAT POS=DEPTH*RAY
   MAT POS=POS+ORIGIN
   IF ABS(DIST)<EPS  THEN EXIT FOR
NEXT I
IF ABS(DIST)<EPS THEN
   CALL GETNORMAL(POS,NORMAL)
   LET DIFFUSE = CLAMP(DOT(LIGHTDIR, NORMAL), 0.1, 1.0 )
   CALL REFLECT(LIGHTDIR, NORMAL,V)
   LET SPECULAR = CLAMP(DOT(V,RAY),0, 1)^10
   CALL VEC3(S,SPECULAR*.8,SPECULAR*.8,SPECULAR*.8)
   MAT P=OFFSET*NORMAL
   MAT P=P+POS
   LET SHADOW = GETSHADOW(P, LIGHTDIR)
   CALL SCENECOLOR(POS ,CL)
   FOR I=1 TO 3
      LET COL(I)=CL(I)
   NEXT I
   MAT COL=DIFFUSE*COL
   MAT COL=COL+S
   MAT COL=(MAX(0.5,SHADOW))*COL
   LET HIT = 1
ELSE
   MAT COL=ZER
   LET HIT=0
END IF
LET K= CLAMP(.05 * DEPTH, 0, .6)^2
FOR I=1 TO 3
   LET COL(I)=COL(I)-K
NEXT I
END SUB

EXTERNAL  SUB REFLECT(I(),N(),V())
DIM C(3),VV(3)
LET K=2*DOT(N,I)
MAT C=K*N
MAT VV=I-C
MAT V=VV
END SUB

EXTERNAL  SUB NORMALIZE(RAY())
LET S=LENGTH(RAY)
IF S<>0 THEN
   MAT RAY=(1/S)*RAY
ELSE
   MAT RAY=ZER
END IF
END SUB

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB VEC2(V(),X,Y)
LET V(1)=X
LET V(2)=Y
END SUB

EXTERNAL  SUB VEC3(V(),X,Y,Z)
LET V(1)=X
LET V(2)=Y
LET V(3)=Z
END SUB

EXTERNAL  SUB VEC4(V(),X,Y,Z,W)
LET V(1)=X
LET V(2)=Y
LET V(3)=Z
LET V(4)=W
END SUB

EXTERNAL  FUNCTION LENGTH(A())
FOR I=1 TO UBOUND(A,1)
   LET S=S+A(I)^2
NEXT I
LET LENGTH=SQR(S)
END FUNCTION

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION FRACT(X)
LET FRACT=FP(X)
END FUNCTION

EXTERNAL  FUNCTION MIX(X,Y,A)
LET MIX=X*(1-A)+Y*A
END FUNCTION



上記プログラムの原版はこちら(※マウスでいじれます)
https://gam0022.net/webgl/#raymarching_reflect

GLSL言語でのGPUによるレンダリングです(WebGL)

※凄すぎるので上記BASICプログラム実行前のアクセス禁止です(笑)
 

レイマーチング

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時13分32秒
返信・引用
  ネット上のサンプルを基にレイマーチング法で鉄骨?をレンダリングしています。
二進モードで実行してください。

DIM P(2),V(2),CPOS(3),CDIR(3),CUP(3),CSIDE(3),RAY(3),R1(3),R2(3),R3(3),LIGHTDIR(3),DPOS(3)
DIM NORMAL(3),COL(3)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
CALL VEC3(CPOS,5,5,.5)
CALL VEC3(CUP,.1,.4,0)
CALL NORMALIZE(CUP)
CALL VEC3(CDIR,-1,0,0)
MAT CDIR=CROSS(CUP,CDIR)
CALL VEC3(LIGHTDIR,1,1,-2)
CALL NORMALIZE(LIGHTDIR)
MAT CSIDE=CROSS(CDIR,CUP)
LET TARGETDEPTH=1
LET EPS=.001
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      MAT R1=X*CSIDE
      MAT R2=Y*CUP
      MAT R3=TARGETDEPTH*CDIR
      MAT RAY=R1
      MAT RAY=RAY+R2
      MAT RAY=RAY+R3
      CALL NORMALIZE(RAY)
      LET DEPTH=0
      MAT DPOS=CPOS
      FOR I=0 TO 63
         LET DIST=DISTANCEFUNC(DPOS)
         LET DEPTH=DEPTH+DIST
         MAT DPOS=DEPTH*RAY
         MAT DPOS=DPOS+CPOS
         IF ABS(DIST)<EPS THEN EXIT FOR
      NEXT I
      IF ABS(DIST)<EPS THEN
         CALL GETNORMAL(DPOS,NORMAL)
         LET DIFFUSE=CLAMP(DOT(LIGHTDIR,NORMAL),.1,1)
         CALL VEC3(COL,1,.1,.1)
         MAT COL=DIFFUSE*COL
         CALL SETCOLOR(COL(1)+.05*DEPTH,COL(2)+.05*DEPTH,COL(3)+.05*DEPTH)
      ELSE
         CALL SETCOLOR(.05*DEPTH,.05*DEPTH,.05*DEPTH)
      END IF
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB  ONREP(P(), INTERVAL,PP())
FOR I=1 TO UBOUND(P,1)
   LET PP(I)=MOD(P(I),INTERVAL)- INTERVAL * 0.5
NEXT I
END SUB

EXTERNAL  FUNCTION BARDIST(X,Y,INTERVAL,WIDTH)
DIM PP(2),P(2)
CALL VEC2(P,X,Y)
CALL ONREP(P, INTERVAL,PP)
CALL VABS(PP)
FOR I=1 TO 2
   LET PP(I)=PP(I)-WIDTH
NEXT I
CALL VMAX(PP,0)
LET BARDIST=LENGTH(PP)
END FUNCTION

EXTERNAL  FUNCTION TUBEDIST(X,Y, INTERVAL, WIDTH)
DIM PP(2),P(2)
CALL VEC2(P,X,Y)
CALL ONREP(P, INTERVAL,PP)
LET TUBEDIST=LENGTH(PP) - WIDTH
END FUNCTION

EXTERNAL  FUNCTION DISTANCEFUNC(P())
LET BARX=BARDIST(P(2),P(3),1,.1)
LET BARY=BARDIST(P(1),P(3),1,.1)
LET BARZ=BARDIST(P(1),P(2),1,.1)
LET TUBEX=TUBEDIST(P(2),P(3),.1,.025)
LET TUBEY=TUBEDIST(P(1),P(3),.1,.025)
LET TUBEZ=TUBEDIST(P(1),P(2),.1,.025)
LET DISTANCEFUNC=MAX(MAX(MAX(MIN(MIN(BARX, BARY),BARZ), -TUBEX), -TUBEY), -TUBEZ)
END FUNCTION

EXTERNAL  SUB GETNORMAL(P(),N())
DIM X1(3),X2(3),Y1(3),Y2(3),Z1(3),Z2(3),PX1(3),PX2(3),PY1(3),PY2(3),PZ1(3),PZ2(3)
LET D=.001
CALL VEC3(X1,D,0,0)
CALL VEC3(X2,-D,0,0)
MAT PX1=P+X1
MAT PX2=P+X2
CALL VEC3(Y1,0,D,0)
CALL VEC3(Y2,0,-D,0)
MAT PY1=P+Y1
MAT PY2=P+Y2
CALL VEC3(Z1,0,0,D)
CALL VEC3(Z2,0,0,-D)
MAT PZ1=P+Z1
MAT PZ2=P+Z2
CALL VEC3(N,DISTANCEFUNC(PX1)-DISTANCEFUNC(PX2),DISTANCEFUNC(PY1)-DISTANCEFUNC(PY2),DISTANCEFUNC(PZ1)-DISTANCEFUNC(PZ2))
CALL NORMALIZE(N)
END SUB

EXTERNAL  SUB NORMALIZE(RAY())
LET S=LENGTH(RAY)
IF S<>0 THEN
   MAT RAY=(1/S)*RAY
ELSE
   MAT RAY=ZER
END IF
END SUB

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB VEC2(V(),X,Y)
LET V(1)=X
LET V(2)=Y
END SUB

EXTERNAL  SUB VEC3(V(),X,Y,Z)
LET V(1)=X
LET V(2)=Y
LET V(3)=Z
END SUB

EXTERNAL  FUNCTION LENGTH(A())
FOR I=1 TO UBOUND(A,1)
   LET S=S+A(I)^2
NEXT I
LET LENGTH=SQR(S)
END FUNCTION

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  SUB VABS(A())
FOR I=1 TO UBOUND(A,1)
   LET A(I)=ABS(A(I))
NEXT I
END SUB

EXTERNAL  SUB VMAX(A(),N)
FOR I=1 TO UBOUND(A,1)
   LET A(I)=MAX(A(I),N)
NEXT I
END SUB


このプログラムの原版はこちら(※マウスでいじれます)
https://gam0022.net/webgl/#raymarching_steel-frame

GLSL言語でのGPUによるレンダリングです(WebGL)

※凄すぎるので上記BASICプログラム実行前のアクセス禁止です(笑)
 

レイマーチング

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時12分42秒
返信・引用
  ネット上のサンプルを基にレイマーチング法によりトーラス(ドーナツ型)をレンダリングしています。
二進モードで実行してください。

DIM CPOS(3),CDIR(3),CUP(3),CSIDE(3),RAY(3),R1(3),R2(3),R3(3),LIGHTDIR(3),RPOS(3),HALF(3)
DIM NORMAL(3),DPOS(3),COL(3)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
CALL VEC3(CPOS,0,5,5)         !'カメラ
CALL VEC3(CDIR,0,-.707,-.707) !'カメラの向き(視線)
CALL VEC3(CUP,0,.707,-.707)   !'カメラの上方向
CALL VEC3(LIGHTDIR,-.577,.577,.577)
MAT CSIDE=CROSS(CDIR,CUP)     !'横方向
LET TARGETDEPTH=1             !'フォーカス深度
CALL NORMALIZE(LIGHTDIR)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      MAT R1=X*CSIDE
      MAT R2=Y*CUP
      MAT R3=TARGETDEPTH*CDIR
      MAT RAY=R1
      MAT RAY=RAY+R2
      MAT RAY=RAY+R3
      CALL NORMALIZE(RAY) !'レイの定義
      LET DIST=0
      LET RLEN=0
      MAT RPOS=CPOS
      LET SHADOW=1
      FOR I=0 TO 255 !'マーチングループ(marching loop)
         CALL DISTANCE(RPOS,DIST,COL)
         LET RLEN=RLEN+DIST
         MAT RPOS=RLEN*RAY
         MAT RPOS=RPOS+CPOS
         IF DIST<.001 THEN EXIT FOR
      NEXT I
      IF ABS(DIST)<.001 THEN !'レイとの距離
         CALL GETNORMAL(RPOS,NORMAL)
         MAT HALF=LIGHTDIR-RAY
         CALL NORMALIZE(HALF)
         LET DIFF=CLAMP(DOT(LIGHTDIR,NORMAL),.1,1)
         LET SPEC=CLAMP(DOT(HALF,NORMAL),0,1)^50
         MAT DPOS=(1/1000)*NORMAL
         MAT DPOS=DPOS+RPOS
         LET SHADOW=GENSHADOW(DPOS,LIGHTDIR) !'シャドウ
         FOR I=1 TO 3
            LET COL(I)=COL(I)*DIFF+SPEC
         NEXT I
      ELSE
         CALL VEC3(COL,0,0,0)
      END IF
      FOR I=1 TO 3
         LET COL(I)=COL(I)*MAX(.5,SHADOW)
      NEXT I
      CALL SETCOLOR(COL(1),COL(2),COL(3))
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  FUNCTION TORUS(P()) !'距離関数(トーラス)
DIM T(3),R(3)
CALL VEC2(T,3,1)
CALL VEC2(R,LENGTH2(P(1),P(3),0)-T(1),P(2))
LET TORUS=LENGTH(R)-T(2)
END FUNCTION

EXTERNAL  FUNCTION FLOOR(P()) !'距離関数(床)
DIM V(3)
CALL VEC3(V,0,1,0)
LET FLOOR=DOT(P,V)+1
END FUNCTION

EXTERNAL  SUB DISTANCE(P(),DIST,COL())
LET D1=TORUS(P)
LET D2=FLOOR(P)
IF D1<D2 THEN
   LET DIST=D1
   CALL VEC3(COL,1,1,.25) !'色
ELSE
   LET DIST=D2
   LET U=1-IP(MOD(P(1),2))
   LET V=1-IP(MOD(P(3),2))
   IF U+V=1 THEN CALL VEC3(COL,.7,.7,.7) ELSE CALL VEC3(COL,1,1,1) !'色(市松模様)
END IF
END SUB

EXTERNAL  SUB NMAX(A(),B(),N)
FOR I=1 TO 3
   LET A(I)=MAX(B(I),N)
NEXT I
END SUB

EXTERNAL  SUB GETNORMAL(P(),N()) !'法線ベクトル
DIM X1(3),X2(3),Y1(3),Y2(3),Z1(3),Z2(3),PX1(3),PX2(3),PY1(3),PY2(3),PZ1(3),PZ2(3),DMY(3)
LET D=.0001
CALL VEC3(X1,D,0,0)
CALL VEC3(X2,-D,0,0)
MAT PX1=P+X1
MAT PX2=P+X2
CALL VEC3(Y1,0,D,0)
CALL VEC3(Y2,0,-D,0)
MAT PY1=P+Y1
MAT PY2=P+Y2
CALL VEC3(Z1,0,0,D)
CALL VEC3(Z2,0,0,-D)
MAT PZ1=P+Z1
MAT PZ2=P+Z2
CALL DISTANCE(PX1,XS,DMY)
CALL DISTANCE(PX2,XE,DMY)
CALL DISTANCE(PY1,YS,DMY)
CALL DISTANCE(PY2,YE,DMY)
CALL DISTANCE(PZ1,ZS,DMY)
CALL DISTANCE(PZ2,ZE,DMY)
CALL VEC3(N,XS-XE,YS-YE,ZS-ZE)
CALL NORMALIZE(N)
END SUB

EXTERNAL  SUB NORMALIZE(RAY()) !'正規化
LET S=LENGTH(RAY)
IF S<>0 THEN
   MAT RAY=(1/S)*RAY
ELSE
   MAT RAY=ZER
END IF
END SUB

EXTERNAL  SUB VABS(A(),B())
FOR I=1 TO 3
   LET A(I)=ABS(B(I))
NEXT I
END SUB

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB VEC3(V(),X,Y,Z)
LET V(1)=X
LET V(2)=Y
LET V(3)=Z
END SUB

EXTERNAL  SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB

EXTERNAL  FUNCTION LENGTH(A()) !'長さ
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION

EXTERNAL  FUNCTION LENGTH2(X,Y,Z)
LET LENGTH2=SQR(X^2+Y^2+Z^2)
END FUNCTION

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION GENSHADOW(RO(),RD()) !'シャドウ
DIM RAY(3),DMY(3)
LET R=1
LET C=.001
LET SHADOWCOEF=.5
FOR T=0 TO 49
   MAT RAY=C*RD
   MAT RAY=RAY+RO
   CALL DISTANCE(RAY,H,DMY)
   IF H<.001 THEN
      LET GENSHADOW=SHADOWCOEF
      EXIT FUNCTION
   END IF
   LET R=MIN(R,H*32/C)
   LET C=C+H
NEXT T
LET GENSHADOW=1-SHADOWCOEF+R*SHADOWCOEF
END FUNCTION
 

レイマーチング

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時11分57秒
返信・引用
  ネット上のサンプルを基にしてレイマーチング法(ray marching)で箱(box)をレンダリングしています。
レイ・マーチング法はレイ・トレーシング法(ray tracing)の一種です。

※厳密にはレイマーチング法の中のスフィアトレーシング法(sphere tracing)です。
二進モードで実行してください。

DIM P(2),V(2),CPOS(3),CDIR(3),CUP(3),CSIDE(3),RAY(3),R1(3),R2(3),R3(3),LIGHTDIR(3),RPOS(3)
DIM NORMAL(3)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
CALL VEC3(CPOS,0,0,2)
CALL VEC3(LIGHTDIR,-.577,.577,.577)
LET ANG=60            !'視野角
LET FOV=ANG*.5*PI/180
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      CALL VEC3(RAY,SIN(FOV)*X,SIN(FOV)*Y,-COS(FOV))
      CALL NORMALIZE(RAY)
      LET DISTANCE=0
      LET RLEN=0
      MAT RPOS=CPOS
      FOR I=0 TO 127
         LET DISTANCE=DISTANCEFUNC(RPOS)
         LET RLEN=RLEN+DISTANCE
         MAT RPOS=RLEN*RAY
         MAT RPOS=RPOS+CPOS
         IF ABS(DISTANCE)<.001 THEN EXIT FOR
      NEXT I
      IF ABS(DISTANCE)<.001 THEN
         CALL GETNORMAL(RPOS,NORMAL)
         LET DIFF=CLAMP(DOT(LIGHTDIR,NORMAL),.1,1)
         CALL SETCOLOR(DIFF,DIFF,DIFF)
      ELSE
         CALL SETCOLOR(0,0,0)
      END IF
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB TRANS(P())
FOR I=1 TO 3
   LET P(I)=MOD(P(I),4)-2
NEXT I
END SUB

EXTERNAL  FUNCTION DISTANCEFUNC(PP())
DIM V(3),A(3),B(3),P(3)
CALL VEC3(B,.5,.5,.5)
MAT P=PP
CALL TRANS(P)
CALL VABS(V,P)
MAT V=V-B
CALL NMAX(A,V,0)
LET DISTANCEFUNC=LENGTH(A)-.5
END FUNCTION

EXTERNAL  SUB GETNORMAL(P(),N())
DIM X1(3),X2(3),Y1(3),Y2(3),Z1(3),Z2(3),PX1(3),PX2(3),PY1(3),PY2(3),PZ1(3),PZ2(3)
LET D=.0001
CALL VEC3(X1,D,0,0)
CALL VEC3(X2,-D,0,0)
MAT PX1=P+X1
MAT PX2=P+X2
CALL VEC3(Y1,0,D,0)
CALL VEC3(Y2,0,-D,0)
MAT PY1=P+Y1
MAT PY2=P+Y2
CALL VEC3(Z1,0,0,D)
CALL VEC3(Z2,0,0,-D)
MAT PZ1=P+Z1
MAT PZ2=P+Z2
CALL VEC3(N,DISTANCEFUNC(PX1)-DISTANCEFUNC(PX2),DISTANCEFUNC(PY1)-DISTANCEFUNC(PY2),DISTANCEFUNC(PZ1)-DISTANCEFUNC(PZ2))
CALL NORMALIZE(N)
END SUB

EXTERNAL  SUB NORMALIZE(RAY())
LET S=LENGTH(RAY)
IF S<>0 THEN
   MAT RAY=(1/S)*RAY
ELSE
   MAT RAY=ZER
END IF
END SUB

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB VEC3(V(),X,Y,Z)
LET V(1)=X
LET V(2)=Y
LET V(3)=Z
END SUB

EXTERNAL  FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2+A(3)^2)
END FUNCTION

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  SUB VABS(A(),B())
FOR I=1 TO 3
   LET A(I)=ABS(B(I))
NEXT I
END SUB

EXTERNAL  SUB NMAX(A(),B(),N)
FOR I=1 TO 3
   LET A(I)=MAX(B(I),N)
NEXT I
END SUB
 

フラグメントシェーダー

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時11分8秒
返信・引用
  フラグメントシェーダーによるプログラムです。
ネット上のサンプルを基にアンパンマン?を描画しています。

https://qiita.com/doxas/items/a366eafc498c8269934c

PUBLIC NUMERIC LIGHTCOLOR(3),BACKCOLOR(3),FACECOLOR(3),NOSECOLOR(3),CHEEKCOLOR(3),EYESCOLOR(3),HIGHLIGHT(3),LINECOLOR(3),T
DIM P(2),M(2,2),V(2),W(2),Q(2),QQ(2),DESTCOLOR(3)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
CALL VEC3(LIGHTCOLOR,0.95, 0.95, 0.5)!'  // 背景の後光の色
CALL VEC3(BACKCOLOR,0.95, 0.25, 0.25)!'  // 背景の下地の色
CALL VEC3(FACECOLOR,0.95, 0.75, 0.5)!'   // 顔の色
CALL VEC3(NOSECOLOR,0.95, 0.25, 0.25)!'  // 鼻の色
CALL VEC3(CHEEKCOLOR,1.0, 0.55, 0.25)!'  // 頬の色
CALL VEC3(EYESCOLOR,0.15, 0.05, 0.05)!'  // 目の色
CALL VEC3(HIGHLIGHT,0.95, 0.95, 0.95)!'  // ハイライトの色
CALL VEC3(LINECOLOR,0.3, 0.2, 0.2)!'     // ラインの色
LET T=INT(TIME)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      CALL VEC2(P,X,Y)
      CALL SUNRISE(P,DESTCOLOR)
      LET S=SIN(SIN(T*2)*.75)
      LET C=COS(SIN(T*2))
      CALL MAT2(M,C,-S,S,C)
      MAT Q=P*M
      !    circle(q, vec2(0.0), 0.5, faceColor, destColor);
      CALL VEC2(V,0,0)
      CALL CIRCLE(Q,V,.5,FACECOLOR,DESTCOLOR)
      !    circle(q, vec2(0.0, -0.05), 0.15, noseColor, destColor);
      CALL VEC2(V,0,-.05)
      CALL CIRCLE(Q,V,.15,NOSECOLOR,DESTCOLOR)
      !    circle(q, vec2(0.325, -0.05), 0.15, cheekColor, destColor);
      CALL VEC2(V,.325,-.05)
      CALL CIRCLE(Q,V,.15,CHEEKCOLOR,DESTCOLOR)
      !    circle(q, vec2(-0.325, -0.05), 0.15, cheekColor, destColor);
      CALL VEC2(V,-.325,-.05)
      CALL CIRCLE(Q,V,.15,CHEEKCOLOR,DESTCOLOR)
      !    ellipse(q, vec2(0.15, 0.135), vec2(0.75, 1.0), 0.075, eyesColor, destColor);
      CALL VEC2(V,.15,.135)
      CALL VEC2(W,.75,1)
      CALL ELLIPSE(Q,V,W,.075,EYESCOLOR,DESTCOLOR)
      !    ellipse(q, vec2(-0.15, 0.135), vec2(0.75, 1.0), 0.075, eyesColor, destColor);
      CALL VEC2(V,-.15,.135)
      CALL ELLIPSE(Q,V,W,.075,EYESCOLOR,DESTCOLOR)
      !    circleLine(q, vec2(0.0), 0.5, 0.525, lineColor, destColor);
      CALL VEC2(V,0,0)
      CALL CIRCLELINE(Q,V,.5,.525,LINECOLOR,DESTCOLOR)
      !    circleLine(q, vec2(0.0, -0.05), 0.15, 0.17, lineColor, destColor);
      CALL VEC2(V,0,-.05)
      CALL CIRCLELINE(Q,V,.15,.17,LINECOLOR,DESTCOLOR)
      !    arcLine(q, vec2(0.325, -0.05), 0.15, 0.17, PI * 1.5, 0.0, lineColor, destColor);
      CALL VEC2(V,.325,-.05)
      CALL ARCLINE(Q,V,.15,.17,PI*.5,0,LINECOLOR,DESTCOLOR)
      !    arcLine(q, vec2(-0.325, -0.05), 0.15, 0.17, PI * 0.5, 0.0, lineColor, destColor);
      CALL VEC2(V,-.325,-.05)
      CALL ARCLINE(Q,V,.15,.17,PI*1.5,0,LINECOLOR,DESTCOLOR)
      !    arcLine(q * vec2(1.2, 1.0), vec2(0.19, 0.2), 0.125, 0.145, 0.0, 0.02, lineColor, destColor);
      CALL VEC2(V,1.2,1)
      CALL VEC2(W,.19,.2)
      FOR I=1 TO 2
         LET QQ(I)=Q(I)*V(I)
      NEXT I
      CALL ARCLINE(QQ,W,.125,.145,0,.02,LINECOLOR,DESTCOLOR)
      !    arcLine(q * vec2(1.2, 1.0), vec2(-0.19, 0.2), 0.125, 0.145, 0.0, 0.02, lineColor, destColor);
      CALL VEC2(V,1.2,1)
      CALL VEC2(W,-.19,.2)
      FOR I=1 TO 2
         LET QQ(I)=Q(I)*V(I)
      NEXT I
      CALL ARCLINE(QQ,W,.125,.145,0,.02,LINECOLOR,DESTCOLOR)
      !    arcLine(q * vec2(0.9, 1.0), vec2(0.0, -0.15), 0.2, 0.22, PI, 0.055, lineColor, destColor);
      CALL VEC2(V,.9,1)
      CALL VEC2(W,0,-.15)
      FOR I=1 TO 2
         LET QQ(I)=Q(I)*V(I)
      NEXT I
      CALL ARCLINE(QQ,W,0.2, 0.22, PI, 0.055,LINECOLOR,DESTCOLOR)
      !    rect(q, vec2(-0.025, 0.0), 0.035, highlight, destColor);
      CALL VEC2(V,-.025,0)
      CALL RECT(Q,V,.035,HIGHLIGHT,DESTCOLOR)
      !    rect(q, vec2(-0.35, 0.0), 0.035, highlight, destColor);
      CALL VEC2(V,-.35,0)
      CALL RECT(Q,V,.035,HIGHLIGHT,DESTCOLOR)
      !    rect(q, vec2(0.3, 0.0), 0.035, highlight, destColor);
      CALL VEC2(V,.3,0)
      CALL RECT(Q,V,.035,HIGHLIGHT,DESTCOLOR)
      CALL SETCOLOR(DESTCOLOR)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB CIRCLE(P(),OFFSET(),SIZE,COL(),I())
DIM PP(2)
MAT PP=P-OFFSET
LET L=LENGTH(PP)
IF L<SIZE THEN MAT I=COL
END SUB

EXTERNAL  SUB ELLIPSE(P(),OFFSET(),PROP(),SIZE,COL(),I())
DIM PP(2),Q(2)
MAT PP=P-OFFSET
FOR J=1 TO 2
   LET Q(J)=PP(J)/PROP(J)
NEXT J
IF LENGTH(Q)<SIZE THEN MAT I=COL
END SUB

EXTERNAL  SUB CIRCLELINE(P(),OFFSET(),ISIZE,OSIZE,COL(),I())
DIM Q(2)
MAT Q=P-OFFSET
LET L=LENGTH(Q)
IF L>ISIZE AND L<OSIZE THEN MAT I=COL
END SUB

EXTERNAL  SUB ARCLINE(P(),OFFSET(),ISIZE,OSIZE,RAD,HEIGHT,COL(),I())
DIM ROT(2,2),Q(2)
LET S=SIN(RAD)
LET C=COS(RAD)
CALL MAT2(ROT,C,-S,S,C)
MAT Q=P-OFFSET
MAT Q=Q*ROT
LET L=LENGTH(Q)
IF L>ISIZE AND L<OSIZE AND Q(2)>HEIGHT THEN MAT I=COL
END SUB

EXTERNAL  SUB RECT(P(),OFFSET(),SIZE,COL(),I())
DIM Q(2)
MAT Q=P-OFFSET
MAT Q=(1/SIZE)*Q
IF ABS(Q(1))<1 AND ABS(Q(2))<1 THEN MAT I=COL
END SUB

EXTERNAL  SUB SUNRISE(P(),I())
LET F=ATAN(P(1),P(2))+T
LET FS=SIN(F*10)
FOR J=1 TO 2
   LET I(J)=MIX(LIGHTCOLOR(J),BACKCOLOR(J),FS)
NEXT J
END SUB

EXTERNAL  SUB SETCOLOR(A())
SET COLOR COLORINDEX(CLAMP(A(1),0,1),CLAMP(A(2),0,1),CLAMP(A(3),0,1))
END SUB

EXTERNAL  SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB

EXTERNAL  SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2)
END FUNCTION

EXTERNAL  FUNCTION MIX(X,Y,A)
LET MIX=X*(1-A)+Y*A
END FUNCTION

EXTERNAL  SUB MAT2(X(,),A,B,C,D)
LET X(1,1)=A
LET X(1,2)=B
LET X(2,1)=C
LET X(2,2)=D
END SUB

EXTERNAL FUNCTION ATAN(X,Y)
IF ABS(X)>1E-4 THEN
   LET TH=ATN(Y/X)
   IF Y<>0 THEN
      IF X>0 AND Y<0 THEN LET TH=TH+PI*2
      IF X<0 THEN LET TH=TH+PI
   ELSE !' Y=0
      IF X<0 THEN LET TH=PI ELSE LET TH=0
   END IF
ELSE !' X=0
   LET TH=PI/2
   IF Y<0 THEN LET TH=TH+PI
END IF
LET ATAN=TH
END FUNCTION
 

フラグメントシェーダー

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時10分27秒
返信・引用
  フラグメントシェーダーによるプログラムです。
ネット上のサンプルを基にきらきら? を描画しています。

https://qiita.com/doxas/items/5d6e39c54e16f352488c

DIM Q(2),M(2,2)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
LET T=INT(TIME)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      LET Q(1)=MOD(X,.2)-.1
      LET Q(2)=MOD(Y,.2)-.1
      LET S=SIN(T)
      LET C=COS(T)
      CALL MAT2(M,C,S,-S,C)
      MAT Q=Q*M
      IF Q(1)*Q(2)=0 THEN
         LET V=1
      ELSE
         LET V=.1/ABS(Q(2))*ABS(Q(1))
      END IF
      LET R=V*ABS(SIN(T*6)+1.5)
      LET G=V*ABS(SIN(T*4.5)+1.5)
      LET B=V*ABS(SIN(T*3)+1.5)
      CALL SETCOLOR(R,G,B)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB MAT2(X(,),A,B,C,D)
LET X(1,1)=A
LET X(1,2)=B
LET X(2,1)=C
LET X(2,2)=D
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION
 

フラグメントシェーダー

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時09分55秒
返信・引用
  フラグメントシェーダーによるプログラムです。
このプログラムでは円を10個描いて花模様?を描いています。

https://qiita.com/doxas/items/25bb50a3db85129e2980

DIM P(2),COL(3),PP(2),Q(2),D(2)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
LET T=INT(TIME)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      CALL VEC2(P,X,Y)
      CALL VEC3(COL,1,.3,.7)
      LET F=0
      FOR I=0 TO 9
         LET S=SIN(T+I*PI/5)*.5
         LET C=COS(T+I*PI/5)*.5
         CALL VEC2(PP,C,S)
         MAT D=P+PP
         LET LL=ABS(LENGTH(D)-.5)
         IF LL<>0 THEN LET F=F+.0025/LL ELSE F=F+1
      NEXT I
      MAT COL=F*COL
      CALL SETCOLOR(COL)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB SETCOLOR(A())
SET COLOR COLORINDEX(CLAMP(A(1),0,1),CLAMP(A(2),0,1),CLAMP(A(3),0,1))
END SUB

EXTERNAL  SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB

EXTERNAL  SUB VEC3(A(),X,Y,Z)
LET A(1)=X
LET A(2)=Y
LET A(3)=Z
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2)
END FUNCTION
 

フラグメントシェーダー

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時09分10秒
返信・引用
  フラグメントシェーダーによるプログラムです。
このプログラムでは十字架?を描画します。

https://qiita.com/doxas/items/25bb50a3db85129e2980

DIM P(2),COL(3),PP(2),Q(2),D(2)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      IF X*Y=0 THEN
         LET F=1
      ELSE
         LET F=.001/ABS(X*Y)
      END IF
      CALL SETCOLOR(F,F,F)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION
 

フラグメントシェーダー

 投稿者:しばっち  投稿日:2018年 4月 9日(月)20時08分37秒
返信・引用
  フラグメントシェーダー(Fragment Shader)では、PLOT LINES文やDRAW CIRCLE文、PLOT AREA文といった描画コマンドは
使用せず、各ピクセル毎に色を定義しながら描画していきます。
その為、描画コマンドは色を決めるSET COLOR命令と点を打つPLOT POINTS文しか使用しません。

ネット上のサンプルを基にGLSL言語(シェーディング言語)を十進BASICに移植してみました。
このプログラムでは光の玉を定義し、描画します。

https://qiita.com/doxas/items/f3f8bf868f12851ea143

DIM P(2)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE) !'範囲を-1~1に正規化(描画領域が正方形の時)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      CALL VEC2(P,X,Y)
      LET L=LENGTH(P)
      IF L<>0 THEN LET C=.1/L ELSE LET C=1
      CALL SETCOLOR(C,C,C)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB SETCOLOR(R,G,B)
SET COLOR COLORINDEX(CLAMP(R,0,1),CLAMP(G,0,1),CLAMP(B,0,1))
END SUB

EXTERNAL  SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2)
END FUNCTION
----------------------------------------------------------------------------------------------------------------
このプログラムでも光の玉(オーブ)を描画しています。
https://qiita.com/doxas/items/00567758621bb506e584

DIM P(2),COL(3),PP(2),Q(2)
ASK BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW -XSIZE/MIN(XSIZE,YSIZE),XSIZE/MIN(XSIZE,YSIZE),-YSIZE/MIN(XSIZE,YSIZE),YSIZE/MIN(XSIZE,YSIZE)
FOR YY=0 TO YSIZE-1
   FOR XX=0 TO XSIZE-1
      MAT COL=ZER
      LET X=(XX*2-XSIZE)/MIN(XSIZE,YSIZE)
      LET Y=(YY*2-YSIZE)/MIN(XSIZE,YSIZE)
      CALL VEC2(P,X,Y)
      FOR I=0 TO 5
         LET J=I+1
         CALL VEC2(PP,COS(J)*.5,SIN(J)*.5)
         MAT Q=P+PP
         LET L=LENGTH(Q)
         FOR K=1 TO 3
            IF L<>0 THEN LET COL(K)=COL(K)+.05/L ELSE LET COL(K)=1
         NEXT K
      NEXT I
      CALL SETCOLOR(COL)
      PLOT POINTS:X,Y
   NEXT XX
NEXT YY
END

EXTERNAL  SUB SETCOLOR(A())
SET COLOR COLORINDEX(CLAMP(A(1),0,1),CLAMP(A(2),0,1),CLAMP(A(3),0,1))
END SUB

EXTERNAL  SUB VEC2(A(),X,Y)
LET A(1)=X
LET A(2)=Y
END SUB

EXTERNAL  FUNCTION CLAMP(X,A,B)
LET CLAMP=MIN(B,MAX(X,A))
END FUNCTION

EXTERNAL  FUNCTION LENGTH(A())
LET LENGTH=SQR(A(1)^2+A(2)^2)
END FUNCTION
 

動作報告です。

 投稿者:たろさ  投稿日:2018年 4月 8日(日)08時58分21秒
返信・引用
  ■ パソコン環境
 ウィンドウズ : Microsoft Windows 10 Pro
 サービスパック : なし
 システムの種類 : 64 ビット
 プロセッサー : Intel(R) Core(TM) i5 CPU       M 560  @ 2.67GHz
 周波数 : 2.67 GHz
 メモリー : 3.73 GB

■ パソコンメーカー
 メーカー : TOSHIBA
 機種名 : dynabook RX3 TN266E/3HD

■ その他
 ハードディスクドライブ 空き容量
  C:ドライブ : 16.1 GB (合計 : 106.3 GB)

 CD/DVD ドライブ
  D:ドライブ : MATSHITA DVD-RAM UJ892ES

 言語設定
  システム言語 : 日本語 (日本)
  ユーザー設定言語 : 日本語 (日本)

 インターネット環境
  Internet Explorer : 11
  ネットワーク接続 : 可能


BASIC783setup.exe (1,577,609 bytes)をインストールさせて頂きました。


行列式、逆行列、連立方程式  投稿者:しばっち  投稿日:2014年11月26日(水)18時10分49秒
http://6317.teacup.com/basic/bbs/3557

掲示板からコピーして実行ボタン ポチットしました。

https://blogs.yahoo.co.jp/donald_stinger

 

レンタル掲示板
/158