久しぶりに更新

夏休みに新しいノートパソコン買って、今まさにcore i5 450Mのスピードに感動しているところだ

最近暇ができてきたのでまたfortranでレイトレしてみた

パソコンも新しくしてスペックも上がり計算に２時間しかかからなかった

これもピクセルごとのRGBをfortranで出力したのをHSPで読み込んでスクリーン出力したものだ

相変わらず進歩しない風景だが一応ソースを載っけてみる

まずはfortranのソース↓（ちなみにコンパイルはg95っての使ってソースのファイル名はh.f90でした）

program IO
integer yoko,tate,kaiso,mate,jimenmate,i,k,m,mcnt1,htmate,mhtmate,hhpi,ii,kk,ll,h,jkl,jjk
real cb0,cb1,cb2,ca0,ca1,ca2,hsa0,hsa1,hsa2,hsb0,hsb1,hsb2,q0,q1,q2,qt2
real t,w,b,c,d4,vmr,vmg,vmb,vmi,kai,x,z,drtga,naiseki,light0,light1,light2
integer mtset(3,640,480)
real g(3,3)
real,allocatable,dimension(:,:) :: kyu
real,allocatable,dimension(:,:) :: ja
real,allocatable,dimension(:,:) :: jb1
real,allocatable,dimension(:,:) :: jb2
real,allocatable,dimension(:,:) :: j12
integer,allocatable,dimension(:,:,:) :: pset

print *,"grnd"
read *, jkl
allocate(ja(3,jkl))
allocate(jb1(3,jkl))
allocate(jb2(3,jkl))
allocate(j12(2,jkl))
jimenmate=jkl

print *,"circle?"
read *, jkl
allocate(kyu(5,jkl))
mate=jkl

print *,"アンチエイリアスx?"
read *, jkl
allocate(pset(3,jkl*640,480*jkl))

print *, "kaiso"
read *, kaiso

yoko=640*jkl
tate=480*jkl
t=-3.0
k=1
w=0.0
jjk=0
do i=1,mate
kyu(4,i)=0.5
kyu(1,i)=w
kyu(2,i)=4.9*mod(i,2)+0.27
kyu(3,i)=t
kyu(5,i)=kyu(4,i)**2
if (mod(i,2).eq.1) then
w=w+0.6
jjk=jjk+1
if (jjk.eq.k) then
jjk=0
w=-0.30*k
k=k+1
t=t+0.99
endif
endif

enddo
w=0.0
t=0.0
jjk=0
k=0

do i=1,jimenmate
ja(1,i)=0.01*rand()*2001-10.0
ja(2,i)=0.01*rand()*601+3.5
ja(3,i)=0.01*rand()*1500+9.25
jb1(1,i)=0.01*rand()*1000-4.95
jb1(2,i)=0.01*rand()*1000-4.95
jb1(3,i)=0.01*rand()*1000-4.95
jb2(1,i)=0.01*rand()*1000-4.95
jb2(2,i)=0.01*rand()*1000-4.95
jb2(3,i)=0.01*rand()*1000-4.95
j12(1,i)=0.7+0.01*rand()*130
j12(2,i)=0.7+0.01*rand()*130
w=sqrt(jb1(1,i)**2.0+jb1(2,i)**2.0+jb1(3,i)**2.0)
jb1(1,i)=jb1(1,i)/w
jb1(2,i)=jb1(2,i)/w
jb1(3,i)=jb1(3,i)/w
w=sqrt(jb2(1,i)**2.0+jb2(2,i)**2.0+jb2(3,i)**2.0)
jb2(1,i)=jb2(1,i)/w
jb2(2,i)=jb2(2,i)/w
jb2(3,i)=jb2(3,i)/w
enddo

light0=1.0*rand()*100-30.0
light1=1.0+rand()*50+rand()*(3.0+rand()*180)
light2=1.0*rand()*100-50.0
w=sqrt(light0**2+light1**2+light2**2)
light0=light0/w
light1=light1/w
light2=light2/w

do k=1,yoko
mcnt1=k
do ii=1,tate
ca0=0.0
ca1=3.0
ca2=-7.0
cb0=-0.004*yoko/jkl+0.008*mcnt1/jkl
cb1=0.004*tate/jkl-0.008*ii/jkl
cb2=9.3
w=sqrt(cb0*cb0+cb1*cb1+cb2*cb2)
cb0=cb0/w
cb1=cb1/w
cb2=cb2/w

hsa0=0.0
hsa1=0.0
hsa2=0.0
hsb0=0.0
hsb1=0.0
hsb2=0.0
mhtmate=-1
vmi=1.0
m=0
do hhpi=1,kaiso
if (m.eq.0) then
htmate=-1
t=99999990.0
do i=1,mate
if (mhtmate.ne.i) then

q0=ca0-kyu(1,i)
q1=ca1-kyu(2,i)
q2=ca2-kyu(3,i)
qt2=q0*q0+q1*q1+q2*q2
b=q0*cb0+q1*cb1+q2*cb2
c=qt2-kyu(5,i)
d4=b*b-c
if (d4.ge.0.0) then
kai=-sqrt(d4)-b
endif
if ((d4.gt.0.0) .and. (kai.gt.0.0) .and. (t.gt.kai)) then
t=kai
htmate=i
endif

endif
enddo
do i=1,jimenmate
if (mhtmate.ne.(i+10000)) then

q0=ca0-ja(1,i)
q1=ca1-ja(2,i)
q2=ca2-ja(3,i)
g(1,1)=cb0
g(2,1)=jb1(1,i)
g(3,1)=jb2(1,i)
g(1,2)=cb1
g(2,2)=jb1(2,i)
g(3,2)=jb2(2,i)
g(1,3)=cb2
g(2,3)=jb1(3,i)
g(3,3)=jb2(3,i)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga.eq.0.0) then
else
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga

if ((kai.lt.j12(1,i)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,i)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
kai=-kai/drtga
if ((kai.gt.0.0).and.(t.gt.kai)) then
t=kai
htmate=i+10000
endif
endif
endif
endif

endif
enddo
if (cb1.lt.0.0) then

x=cb0*ca1/abs(cb1)
z=cb2*ca1/abs(cb1)
kai=sqrt(x*x+ca1*ca1+z*z)
x=x+ca0
z=z+ca2

if (kai.lt.t) then
if (mod(abs(int(x*2.0)+int(z*2.0)),2) .eq. 1) then
vmr=10.0
vmg=10.0
vmb=93.0
else
vmr=155.0
vmg=123.0
vmb=190.0
endif
w=abs(x)+abs(z)
if (w.lt.400.0) then
vmr=vmr*(400.0-w)/400.0+3.0
vmb=vmb*(400.0-w)/400.0+4.0
vmg=vmg*(400.0-w)/400.0+11.0
else
vmr=3.0
vmb=4.0
vmg=11.0
endif
ca0=x
ca1=0.0
ca2=z
x=1.0
do kk=1,mate
q0=ca0-kyu(1,kk)
q1=-kyu(2,kk)
q2=ca2-kyu(3,kk)
qt2=q0*q0+q1*q1+q2*q2
b=q0*light0+q1*light1+q2*light2
c=qt2-kyu(5,kk)
if ((b*b-c).gt.(0.0)) then
x=0.0
endif
enddo

do kk=1,jimenmate
q0=ca0-ja(1,kk)
q1=-ja(2,kk)
q2=ca2-ja(3,kk)
g(1,1)=light0
g(2,1)=jb1(1,kk)
g(3,1)=jb2(1,kk)
g(1,2)=light1
g(2,2)=jb1(2,kk)
g(3,2)=jb2(2,kk)
g(1,3)=light2
g(2,3)=jb1(3,kk)
g(3,3)=jb2(3,kk)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga .ne. 0.0) then
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga
if ((kai.lt.j12(1,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
if ((kai/drtga) .lt. 0.0) then
x=0.0
endif
endif
endif
endif
enddo

if (x.eq.1.0) then
naiseki=cb0*light0-cb1*light1+cb2*light2
vmr=vmr+(1.3+naiseki)*55.0
if (naiseki.gt.0.95) then
vmr=vmr+(naiseki-0.95)*1730.0
endif
vmg=vmg+(1.3+naiseki)*35.0
if (naiseki.gt.0.95) then
vmg=vmg+(naiseki-0.95)*900.0
endif
vmb=vmb+(1.3+naiseki)*17.16
if (naiseki.gt.0.95) then
vmb=vmb+(naiseki-0.95)*450.0
endif
else
vmr=vmr/1.7
vmg=vmr/1.7
vmb=vmr/1.7
endif
htmate=-1
endif

endif
mhtmate=htmate
if (htmate.eq.-1) then
m=1
else

if (htmate.lt.10000) then
hsb0=(ca0-kyu(1,htmate)+t*cb0)/kyu(4,htmate)
hsb1=(ca1-kyu(2,htmate)+t*cb1)/kyu(4,htmate)
hsb2=(ca2-kyu(3,htmate)+t*cb2)/kyu(4,htmate)
hsa0=kyu(1,htmate)+hsb0*kyu(4,htmate)
hsa1=kyu(2,htmate)+hsb1*kyu(4,htmate)
hsa2=kyu(3,htmate)+hsb2*kyu(4,htmate)
endif

if (htmate.ge.10000) then
htmate=htmate-10000
hsa0=ca0+t*cb0
hsa1=ca1+t*cb1
hsa2=ca2+t*cb2
hsb0=jb1(2,htmate)*jb2(3,htmate)-jb1(3,htmate)*jb2(2,htmate)
hsb1=jb1(3,htmate)*jb2(1,htmate)-jb1(1,htmate)*jb2(3,htmate)
hsb2=jb1(1,htmate)*jb2(2,htmate)-jb1(2,htmate)*jb2(1,htmate)
w=sqrt(hsb0*hsb0+hsb1*hsb1+hsb2*hsb2)
hsb0=hsb0/w
hsb1=hsb1/w
hsb2=hsb2/w
htmate=htmate+10000
endif

naiseki=-(hsb0*cb0+hsb1*cb1+hsb2*cb2)
ca0=hsa0
ca1=hsa1
ca2=hsa2
cb0=2.0*naiseki*hsb0+cb0
cb1=2.0*naiseki*hsb1+cb1
cb2=2.0*naiseki*hsb2+cb2
vmi=vmi+0.7

endif
endif

enddo

if (cb1.ge.0.0) then

naiseki=(cb0*light0+cb1*light1+cb2*light2)
vmr=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmr=vmr+(naiseki-0.96)*8200.0
endif
vmg=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmg=vmg+(naiseki-0.96)*8200.0
endif
vmb=(1.5+naiseki)*190.0
if (naiseki.gt.0.96) then
vmb=vmb+(naiseki-0.96)*8200.0
endif

endif

ll=int(vmr/vmi)
if (ll.lt.256) then
pset(3,mcnt1,ii)=ll
else
pset(3,mcnt1,ii)=255
endif

ll=int(vmg/vmi)
if (ll.lt.256) then
pset(2,mcnt1,ii)=ll
else
pset(2,mcnt1,ii)=255
endif

ll=int(vmb/vmi)
if (ll.lt.256) then
pset(1,mcnt1,ii)=ll
else
pset(1,mcnt1,ii)=255
endif

vmr=0.0
vmg=0.0
vmb=0.0
enddo
enddo

open(11,file = 'file1.txt')
open(12,file = 'file2.txt')
open(13,file = 'file3.txt')
open(14,file = 'file4.txt')
i=0
h=0
k=tate/jkl+1
m=0
ii=0
kk=1
ll=0
do i=1,yoko
do h=1,tate
do m=1,3
mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)=mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)+pset(m,i,h)
enddo
enddo
enddo
m=0
h=0
i=0

hhpi=yoko*tate/jkl/jkl*3/4
do i=1,yoko*tate/jkl/jkl*3
h=mod(i,3)+1
m=mod(i,4)+1
if (h.eq.1) then
ii=mod(ii,yoko/jkl)+1
if (ii.eq.1) then
k=k-1
endif
endif
ll=ll+(mtset(h,ii,k)/jkl/jkl)*kk
kk=kk*256
if (m.eq.4) then
write(11+(i-1)/hhpi,*) ll
ll=0
kk=1
endif
enddo
close(11)
close(12)
close(13)
close(14)
print *,"end111111111111111111111"
stop
end

私は実行するときコマンドプロンプトから実行しているが、多分これを実行すると最初に

grnd?

と入力を求められる

ここに「四角形の板」の枚数を入力する。鏡みたいなやつだ。当然反射する。

上の生成画像では確か１０って設定だったかな？

次に

circle?

と入力を求められる

ここに「球」の個数を入力する

上の画像では確か15000くらい。

アンチエイリアス?

にはアンチエイリアスなしで1を、2,3,4となるごとにきめ細かくなっていくが計算時間が２乗倍になる

上の画像では確か4

kaiso?

は、最高反射計算回数。まぁ4以上の数字ならたいてい変な画像はできない

上の画像は確か19とか

これでenterキーを押せば計算開始、計算が終了する頃に「end111111111111111111111」っていう計算終了のお知らせメッセージが表示される

また出力ファイルで４つのテキストファイルが出力される

そこに全ピクセルのRGB情報が詰まっている

今度はHSPでスクリーン出力

そのソースは↓

onexit*enf
repeat 4
ww=3-cnt
exist ""+ww+""
if strsize=-1:IIIIIIIIIIDDDDDDD=ww:bsave ""+ww+"",IIIIIIIIIIDDDDDDD:break
if cnt=3:end
loop

screen 0,640,480:boxf
dim i,640*480*3/4
mref i,66
sdim k,16

h=640*480*3/4/4*IIIIIIIIIIDDDDDDD
repeat 1,IIIIIIIIIIDDDDDDD
sdim a,1
notesel a
noteload "file"+(cnt+1)+".txt"
repeat 640*480*3/4/4
noteget k,cnt
i.h=int(k)
h+
if cnt\900=0:await 1:redraw
loop
loop

bmpsave ""+IIIIIIIIIIDDDDDDD+".bmp"
delete ""+IIIIIIIIIIDDDDDDD+""
if IIIIIIIIIIDDDDDDD=1{
repeat 999999
w=0
repeat 4
exist ""+cnt+""
if strsize!-1:w+
loop
if w=0:break
wait 200
loop
screen 0,640,480:boxf
repeat 4
buffer 1+cnt,640,480:picload ""+cnt+".bmp"
loop
gsel 0
repeat 4
pos 0,0
gmode 5,640,480,256
pos 0,0
gcopy 1+cnt,0,0,640,480
loop
bmpsave "5.bmp"
delete "0.bmp"
delete "1.bmp"
delete "2.bmp"
delete "3.bmp"
stop
}else{
end
}
*enf
delete ""+IIIIIIIIIIDDDDDDD+""
end

相変わらず超汚い

こんなの誰も理解できんだろうな・・そのうち自分も分からなくなるよこれ

実はこのプログラムDualコア4スレッド対応プログラム

でもHSPはシングルスレッドなのでやっぱ４回起動する必要があるｗ

１スレッドで普通にプログラム実行すると1分くらいかかるので折角のcore i5なので「並列プログラム」を組んでみたというわけ

注意点はなるべく素早く４回起動しないとバグる！(ぇ

しかし別に原理的にはただ単に４つ起動させているだけなので１コアのCPUでも全然動くｗ

で、４つ起動させて全部の計算が終わると残り３つのプログラムは自動終了して１つのスクリーンに完成画像が表示されて「5.bmp」ってファイルに出力される。

これでやっと上の画像が完成

おしまい

っと、今のところfortranとHSPの組み合わせ研究はこんな感じで進んでいる。

今度はfortranでムービーの作成でもやってみようかと考えている！

でRGBデータ→ムービーの変換のところでHSPを使うって感じで