リンクその２

リンクその２：影のある積み木プロジェクト

[コード](2012/12/8改訂)

! Neo cubes ver 6.39
! (C) 1985 - 2012 by （名前を入れてください）
! 横表示用
! オプション /Qip /O3 /Qparallel /QxHost
! リンカー追加ライブラリー mkl_solver.lib mkl_intel_c_dll.lib mkl_intel_thread_dll.lib mkl_core_dll.lib mkl_lapack95.lib
! 「面のグループ化」アルゴリズム実験（建物に有利）maxcapture 要調整
! 長い影を切り取り正しく表示する
! 表示面を切り取り正しく表示する

include 'lapack.f90' !LAPACKを使用する

!共通変数
module common_variables

use ifqwin
use lapack95 !LAPACKを使用する
use f95_precision !LAPACKを使用する

!Graphic画面の大きさを定数で定義する(1600 x 900)
integer(2), parameter :: size_x=1600, size_y=900

!チャイルドウインドウの装置番号
integer(4), parameter :: defaultwindow=0, childwindow1=2

!テキスト出力ウインドウの装置番号
integer(4), parameter :: textwindow=2
!テキスト出力横座標
integer(4), parameter :: text_x=1

!各種受け渡し変数
integer(4) result
integer(2) result2
type (xycoord) xy
TYPE (wxycoord) wxy
TYPE (rccoord) curpos
integer(2) key

!仮想画面 (x, y)
integer(4) vram(0:size_x-1, 0:size_y-1)
integer(4) vram2(0:size_x-1, 0:size_y-1)
integer(4) vram3(0:size_x-1, 0:size_y-1)

!四角リージョン変数 (0, 0)-(size_x-1, size_y-1) のリージョンが最初に１度だけ設定される
real(8) rctx(0:3), rcty(0:3) !頂点の座標（定数）
real(8) rctrgn_a(0:3), rctrgn_b(0:3), rctrgn_c(0:3) !各辺のabc
real(8) rctrgn_sign(0:3) !ポリゴンリージョンの真ん中の点の各辺に対する符号

!３角形表示作業用変数
integer(2) trix(0:2), triy(0:2)

!リージョンに表示しようとするポリゴンの頂点
integer(2) dots !poly()用
real(8) plx(0:31), ply(0:31)
type (xycoord) dispdata(0:31)
integer(2) rgnxm, rgnym, rgnxp, rgnyp
integer(2) rgnxm2, rgnym2, rgnxp2, rgnyp2
integer(2) afterpolydots

!交点計算用入力変数
real(8) crx(0:3), cry(0:3)

!頂点、交点のバッファ
real(8) cxbuf(0:31), cybuf(0:31)
integer(2) ptrcbuf(0:31)
real(8) degree_cbuf(0:31)
integer(2) spcbuf

!getabc用変数
real(8) line_x0, line_y0, line_x1, line_y1
real(8) line_a, line_b, line_c

!ポリゴン表示用変数（３～３２角形）
TYPE (wxycoord) polydata(0:31)

!getabcd用変数
real(8) sur_x0, sur_y0, sur_z0, sur_x1, sur_y1, sur_z1, sur_x2, sur_y2, sur_z2
real(8) sur_a, sur_b, sur_c, sur_d

!立方体の色
integer(4) cubecolor(0:5)
!影の色(and処理)
integer(4), parameter :: shadowcolor=16#7f7f7f

!回転前の面オフセット（定数） cubeoffset(xyz, 頂点番号, 面番号)
real(8) cubeoffset(0:2, 0:3, 0:5)
!回転後の面オフセット rotoffset(xyz, 頂点番号, 面番号)
real(8) rotoffset(0:2, 0:3, 0:5)
!接する面の方向オフセット（定数） dircontoffset(xyz, 面番号)
integer(4) dircontoffset(0:2, 0:5)
!回転用グローバル変数
real(8) rxx, ryy, rzz

!視点の位置
integer(4) eye(0:2)

!視点の方向
real(8) matrix0(0:2, 0:2)
real(8) matrix(0:2, 0:2)
real(8) seev(0:2)

!回転用定数
real(8) sin_t, cos_t

!光源のベクトル
real(8) ray(0:2)
real(8) rotray(0:2) !回転後の光源ベクトル

!立方体の座標データ
integer(4), parameter :: maxcubes=10000
integer(4) cube(0:2, 0:maxcubes-1) !cube(xyz, 通し番号)
integer(4) ptrcube !現在の立方体の個数
real(8) rot(0:2, 0:maxcubes-1) !回転後のデータ
real(8) rot_a(0:5), rot_b(0:5), rot_c(0:5) !abc(面番号)
real(8) rot_d(0:5, 0:maxcubes-1) !d(面番号, 通し番号)
logical(4) act(0:5, 0:maxcubes-1) !act(面番号, 通し番号)
logical(4) actshadow(0:5) !actshadow(面番号)
logical(4) actcontact(0:5, 0:maxcubes-1) !actcontact(面番号, 通し番号)互いに接している面は除外する
real(8) sur(0:2, 0:3, 0:5, 0:maxcubes-1) !sur(xyz, 頂点番号, 面番号, 通し番号)頂点の座標
integer(4) sortptr(0:maxcubes-1)
integer(4) endofsortptr !現在のソートポインタの個数
real(8) distance(0:maxcubes-1)
real(8) signc(0:5) !signc(面番号)立方体の真ん中の点について、６つの面に対する正負を計算
real(8) centerpoint(0:2, 0:5, 0:maxcubes-1) !centerpoint(xyz, 面番号, 通し番号) 各面の真ん中の点

!面のグループ化
integer(4) ppgroup, pgroupforcube
real(8) groupd(0:maxcubes*3-1)
integer(4) groupforsurface(0:maxcubes*3-1)
integer(4) ptrgroup(0:maxcubes*3-1)
integer(4) groupforserial(0:maxcubes*3-1)
integer(4) endofnext(0:maxcubes*3-1)
integer(4) nextptr(0:maxcubes*3-1)
integer(4) ppsort
integer(4) sortc(0:maxcubes*3-1)
logical(4) identity(0:maxcubes-1)
logical(4) isrgnsetted
integer(4), parameter :: maxcapture=maxcubes*8000 !*******要調整
integer(4) rgnxmgroup, rgnxpgroup, rgnymgroup, rgnypgroup
integer(4) capture(0:maxcapture-1)

!カットライン付近の処理
integer(4) ptrcut(0:5, 0:maxcubes-1) !ptrcut(面番号, 通し番号) 各立方体の各面について、カットライン付近の切り取り処理後の座標データがある場合はそのポインタを返す。無ければ-1
integer(4) endofcutangles !cutanglesへのポインタ
integer(4), parameter :: maxcutangles=20000 !*******要調整
real(8) cutangles(0:maxcutangles-1) !頂点の数、頂点x0,y0, x1,y1, ・・・のデータ
real(8), parameter :: sqr3d2=0.867_8 !sqr(3)/2

!影の表示モード
logical(4) shadowmode
logical(4) detailmode

!スクリーンのZ座標
real(8) cutline, screen, eyeoffset
real(8) cutline2 !Detail mode用
integer(1) screenmode

!表示倍率
integer(4) zoom

!回転１ストロークあたりの分割数（大きくすると、なめらかに回転する）
integer(4), parameter :: rollingdiv=10

!ファイル操作
!ファイル読み出し、書き込み用装置番号
integer(4), parameter :: fileequipnumber=11
integer(2) filenumber !ファイル番号
character(10) filenamelist(35)

!表示面のカウンタ（デバッグ用）*******
integer(4) surfacecount, shadowcount

logical(4) debugmode
integer(2) capturepercent
integer(2) cutanglespercent

logical(4) helpmode

end module common_variables

!メインプログラム
use common_variables
implicit none

integer(4) j, k, co
integer(1) d1,d2,d3
integer(4) d
integer(4) the_offset, the_ptr

real(8) clock_start, clock_stop, dclock

INTEGER(1), ALLOCATABLE:: gbuf (:)
INTEGER(2) error
INTEGER(4) imsize

imsize = IMAGESIZE (INT2(0), INT2(0), INT2(size_x-1), INT2(size_y-1))
ALLOCATE(gbuf(imsize), STAT = error)
IF (error .NE. 0) THEN
STOP 'ERROR: Insufficient memory'
END IF

!デバッグモードを初期化する
debugmode=.false.

!ヘルプ表示モードを初期化する
helpmode=.false.

!vram2, vram3 をクリアする
vram2=0
vram3=16#ffffff

!四角リージョンを初期設定する
call setrctrgn()

!ファイル名を初期設定する
filenamelist(1)='cube1.txt'
filenamelist(2)='cube2.txt'
filenamelist(3)='cube3.txt'
filenamelist(4)='cube4.txt'
filenamelist(5)='cube5.txt'
filenamelist(6)='cube6.txt'
filenamelist(7)='cube7.txt'
filenamelist(8)='cube8.txt'
filenamelist(9)='cube9.txt'
filenamelist(10)='cube10.txt' !a
filenamelist(11)='cube11.txt' !b
filenamelist(12)='cube12.txt' !c
filenamelist(13)='cube13.txt' !d
filenamelist(14)='cube14.txt' !e
filenamelist(15)='cube15.txt' !f
filenamelist(16)='cube16.txt' !g
filenamelist(17)='cube17.txt' !h
filenamelist(18)='cube18.txt' !i
filenamelist(19)='cube19.txt' !j
filenamelist(20)='cube20.txt' !k
filenamelist(21)='cube21.txt' !l
filenamelist(22)='cube22.txt' !m
filenamelist(23)='cube23.txt' !n
filenamelist(24)='cube24.txt' !o
filenamelist(25)='cube25.txt' !p
filenamelist(26)='cube26.txt' !q
filenamelist(27)='cube27.txt' !r
filenamelist(28)='cube28.txt' !s
filenamelist(29)='cube29.txt' !t
filenamelist(30)='cube30.txt' !u
filenamelist(31)='cube31.txt' !v
filenamelist(32)='cube32.txt' !w
filenamelist(33)='cube33.txt' !x
filenamelist(34)='cube34.txt' !y
filenamelist(35)='cube35.txt' !z

!影表示モードを初期設定する
shadowmode=.false.
detailmode=.false.

!立方体の個数＝０と初期設定する
ptrcube=0

!スクリーンモード設定
screenmode=1 !1 または 2
call setcutline()

!sin_t, cos_t（定数）を設定する
sin_t= dsin(datan(1.0_8)/45*90/rollingdiv)
cos_t= dcos(datan(1.0_8)/45*90/rollingdiv)

!視点を初期設定する
!位置
eye(0)=10000
eye(1)=10000
eye(2)=10000
!方向
call reseteye001()

!光源ベクトルを初期設定する
ray(0)=0.0_8
ray(1)=0.0_8
ray(2)=1.0_8

!オフセット（定数）を設定する
!面番号0 = x-
cubeoffset(0,0,0)=-0.5_8 !x
cubeoffset(1,0,0)=-0.5_8 !y
cubeoffset(2,0,0)=-0.5_8 !z
cubeoffset(0,1,0)=-0.5_8 !x
cubeoffset(1,1,0)=+0.5_8 !y
cubeoffset(2,1,0)=-0.5_8 !z
cubeoffset(0,2,0)=-0.5_8 !x
cubeoffset(1,2,0)=+0.5_8 !y
cubeoffset(2,2,0)=+0.5_8 !z
cubeoffset(0,3,0)=-0.5_8 !x
cubeoffset(1,3,0)=-0.5_8 !y
cubeoffset(2,3,0)=+0.5_8 !z
dircontoffset(0,0)=-1
dircontoffset(1,0)=0
dircontoffset(2,0)=0
!面番号1 = x+
cubeoffset(0,0,1)=+0.5_8 !x
cubeoffset(1,0,1)=-0.5_8 !y
cubeoffset(2,0,1)=-0.5_8 !z
cubeoffset(0,1,1)=+0.5_8 !x
cubeoffset(1,1,1)=+0.5_8 !y
cubeoffset(2,1,1)=-0.5_8 !z
cubeoffset(0,2,1)=+0.5_8 !x
cubeoffset(1,2,1)=+0.5_8 !y
cubeoffset(2,2,1)=+0.5_8 !z
cubeoffset(0,3,1)=+0.5_8 !x
cubeoffset(1,3,1)=-0.5_8 !y
cubeoffset(2,3,1)=+0.5_8 !z
dircontoffset(0,1)=+1
dircontoffset(1,1)=0
dircontoffset(2,1)=0
!面番号2 = y-
cubeoffset(0,0,2)=-0.5_8 !x
cubeoffset(1,0,2)=-0.5_8 !y
cubeoffset(2,0,2)=-0.5_8 !z
cubeoffset(0,1,2)=+0.5_8 !x
cubeoffset(1,1,2)=-0.5_8 !y
cubeoffset(2,1,2)=-0.5_8 !z
cubeoffset(0,2,2)=+0.5_8 !x
cubeoffset(1,2,2)=-0.5_8 !y
cubeoffset(2,2,2)=+0.5_8 !z
cubeoffset(0,3,2)=-0.5_8 !x
cubeoffset(1,3,2)=-0.5_8 !y
cubeoffset(2,3,2)=+0.5_8 !z
dircontoffset(0,2)=0
dircontoffset(1,2)=-1
dircontoffset(2,2)=0
!面番号3 = y+
cubeoffset(0,0,3)=-0.5_8 !x
cubeoffset(1,0,3)=+0.5_8 !y
cubeoffset(2,0,3)=-0.5_8 !z
cubeoffset(0,1,3)=+0.5_8 !x
cubeoffset(1,1,3)=+0.5_8 !y
cubeoffset(2,1,3)=-0.5_8 !z
cubeoffset(0,2,3)=+0.5_8 !x
cubeoffset(1,2,3)=+0.5_8 !y
cubeoffset(2,2,3)=+0.5_8 !z
cubeoffset(0,3,3)=-0.5_8 !x
cubeoffset(1,3,3)=+0.5_8 !y
cubeoffset(2,3,3)=+0.5_8 !z
dircontoffset(0,3)=0
dircontoffset(1,3)=+1
dircontoffset(2,3)=0
!面番号4 = z-
cubeoffset(0,0,4)=-0.5_8 !x
cubeoffset(1,0,4)=-0.5_8 !y
cubeoffset(2,0,4)=-0.5_8 !z
cubeoffset(0,1,4)=+0.5_8 !x
cubeoffset(1,1,4)=-0.5_8 !y
cubeoffset(2,1,4)=-0.5_8 !z
cubeoffset(0,2,4)=+0.5_8 !x
cubeoffset(1,2,4)=+0.5_8 !y
cubeoffset(2,2,4)=-0.5_8 !z
cubeoffset(0,3,4)=-0.5_8 !x
cubeoffset(1,3,4)=+0.5_8 !y
cubeoffset(2,3,4)=-0.5_8 !z
dircontoffset(0,4)=0
dircontoffset(1,4)=0
dircontoffset(2,4)=-1
!面番号5 = z+
cubeoffset(0,0,5)=-0.5_8 !x
cubeoffset(1,0,5)=-0.5_8 !y
cubeoffset(2,0,5)=+0.5_8 !z
cubeoffset(0,1,5)=+0.5_8 !x
cubeoffset(1,1,5)=-0.5_8 !y
cubeoffset(2,1,5)=+0.5_8 !z
cubeoffset(0,2,5)=+0.5_8 !x
cubeoffset(1,2,5)=+0.5_8 !y
cubeoffset(2,2,5)=+0.5_8 !z
cubeoffset(0,3,5)=-0.5_8 !x
cubeoffset(1,3,5)=+0.5_8 !y
cubeoffset(2,3,5)=+0.5_8 !z
dircontoffset(0,5)=0
dircontoffset(1,5)=0
dircontoffset(2,5)=+1

!色を設定する
co=255-32*0 !x-
cubecolor(0)=co+256*co+65536*co
co=255-32*0 !x+
cubecolor(1)=co+256*co+65536*co
co=255-32*1 !y-
cubecolor(2)=co+256*co+65536*co
co=255-32*1 !y+
cubecolor(3)=co+256*co+65536*co
co=255-32*2 !z-
cubecolor(4)=co+256*co+65536*co
co=255-32*2 !z+
cubecolor(5)=co+256*co+65536*co

OPEN(childwindow1, file='user', title='Graphics child window-1')

!チャイルドウインドウの大きさを設定する
call setchildwindowwidth()

result = SETACTIVEQQ (defaultwindow)

!画面をクリアする
result=setcolorrgb(16#000000)
result2=RECTANGLE( $GFILLINTERIOR, int2(0), int2(0), int2(size_x-1), int2(size_y-1) )
!クリア画面をそのまま取り込む（最初の謎の５４バイトを取得）
CALL GETIMAGE (INT2(0), INT2(0), INT2(size_x-1), INT2(size_y-1), gbuf)

111 continue

key=16#0000

do while ( (key .ne. ichar('q')) .and. (key .ne. ichar('Q')) )

if(debugmode) clock_start=dclock()

!vramをクリアする
vram=0

!表示する
if(shadowmode .and. detailmode) then
call displaydetail()
else
call displaycube()
end if

!データを変換する
the_offset=54
the_ptr=the_offset
do j=0, size_y-1
do k=0, size_x-1
d=vram(k,j)
d1= (d/65536 .and. 255)
d2= (d/256 .and. 255)
d3= (d .and. 255)
gbuf(the_ptr+0)=0
gbuf(the_ptr+1)=d1
gbuf(the_ptr+2)=d2
gbuf(the_ptr+3)=d3
the_ptr= the_ptr+4
end do
end do

!転送する
result = SETACTIVEQQ (childwindow1)
call putimage( int2(0), int2(0), gbuf, $GPSET )

if(debugmode) clock_stop=dclock()

!フォーカス（見せるウインドウ）をチャイルドウインドウ１に設定
result=focusqq(childwindow1)

result = SETACTIVEQQ (childwindow1)

if(helpmode) then !ヘルプを表示

!文字の出力位置を設定する（Y座標、X座標の順）
CALL SETTEXTPOSITION (INT2(1), INT2(text_x), curpos)
write(childwindow1,*) '<-------'
CALL SETTEXTPOSITION (INT2(2), INT2(text_x), curpos)
write(childwindow1,*) ' Neo cubes shader'
CALL SETTEXTPOSITION (INT2(3), INT2(text_x), curpos)
write(childwindow1,*) ' ------->'
CALL SETTEXTPOSITION (INT2(5), INT2(text_x), curpos)
write(childwindow1,*) 'key function:'
CALL SETTEXTPOSITION (INT2(6), INT2(text_x), curpos)
write(childwindow1,*) ' 2,4,6,8,x,y=rotate'
CALL SETTEXTPOSITION (INT2(7), INT2(text_x), curpos)
write(childwindow1,*) ' z,a,s,w,spc,b=move'
CALL SETTEXTPOSITION (INT2(8), INT2(text_x), curpos)
write(childwindow1,*) ' 1=create/delete'
CALL SETTEXTPOSITION (INT2(9), INT2(text_x), curpos)
write(childwindow1,*) ' m=shadow mode'
CALL SETTEXTPOSITION (INT2(10), INT2(text_x), curpos)
write(childwindow1,*) ' r=ray setting'
CALL SETTEXTPOSITION (INT2(11), INT2(text_x), curpos)
write(childwindow1,*) ' i=load file'
CALL SETTEXTPOSITION (INT2(12), INT2(text_x), curpos)
write(childwindow1,*) ' o=save file'
CALL SETTEXTPOSITION (INT2(13), INT2(text_x), curpos)
write(childwindow1,*) ' e=reset eye'
CALL SETTEXTPOSITION (INT2(14), INT2(text_x), curpos)
write(childwindow1,*) ' c=screen mode'
CALL SETTEXTPOSITION (INT2(15), INT2(text_x), curpos)
write(childwindow1,*) ' /=detail mode'

!文字ワードラップをoffにする
result2 = WRAPON($GWRAPOFF)

CALL SETTEXTPOSITION (INT2(16), INT2(text_x), curpos)
write(childwindow1,*) 'cubes:', ptrcube
CALL SETTEXTPOSITION (INT2(17), INT2(text_x), curpos)
write(childwindow1,*) 'pos x:', eye(0)
CALL SETTEXTPOSITION (INT2(18), INT2(text_x), curpos)
write(childwindow1,*) 'pos y:', eye(1)
CALL SETTEXTPOSITION (INT2(19), INT2(text_x), curpos)
write(childwindow1,*) 'pos z:', eye(2)
CALL SETTEXTPOSITION (INT2(20), INT2(text_x), curpos)
write(childwindow1,*) 'surfaces:', surfacecount
CALL SETTEXTPOSITION (INT2(21), INT2(text_x), curpos)
write(childwindow1,*) 'eye direction:'
CALL SETTEXTPOSITION (INT2(22), INT2(text_x), curpos)
write(childwindow1,*) real(matrix0(2,0))
CALL SETTEXTPOSITION (INT2(23), INT2(text_x), curpos)
write(childwindow1,*) real(matrix0(2,1))
CALL SETTEXTPOSITION (INT2(24), INT2(text_x), curpos)
write(childwindow1,*) real(matrix0(2,2))

if(shadowmode) then
CALL SETTEXTPOSITION (INT2(25), INT2(text_x), curpos)
write(childwindow1,*) 'ray direction:'
CALL SETTEXTPOSITION (INT2(26), INT2(text_x), curpos)
write(childwindow1,*) real(ray(0))
CALL SETTEXTPOSITION (INT2(27), INT2(text_x), curpos)
write(childwindow1,*) real(ray(1))
CALL SETTEXTPOSITION (INT2(28), INT2(text_x), curpos)
write(childwindow1,*) real(ray(2))
CALL SETTEXTPOSITION (INT2(29), INT2(text_x), curpos)
write(childwindow1,*) 'shadows:', shadowcount
if(detailmode) then
CALL SETTEXTPOSITION (INT2(30), INT2(text_x), curpos)
write(childwindow1,*) 'Detail mode'
end if
end if

if(debugmode) then
CALL SETTEXTPOSITION (INT2(31), INT2(text_x), curpos)
write(childwindow1,*) 'time=', real(clock_stop-clock_start)
end if

if( debugmode .and. shadowmode .and. (.not.detailmode) ) then
CALL SETTEXTPOSITION (INT2(32), INT2(text_x), curpos)
write(childwindow1,*) 'capture=', capturepercent, '%'
CALL SETTEXTPOSITION (INT2(33), INT2(text_x), curpos)
write(childwindow1,*) 'cutline=', cutanglespercent, '%'
end if

CALL SETTEXTPOSITION (INT2(34), INT2(text_x), curpos)
write(childwindow1,*) 'screen mode=', screenmode

else

!文字の出力位置を設定する（Y座標、X座標の順）
CALL SETTEXTPOSITION (INT2(1), INT2(text_x), curpos)
write(childwindow1,*) 'h=help'

end if

!フォーカス（見せるウインドウ）をチャイルドウインドウ１に設定
result=focusqq(childwindow1)

key=incharqq()
select case (key)
case(ichar('4'))
call rolling(2, 0, +1)
case(ichar('6'))
call rolling(2, 0, -1)
case(ichar('8'))
call rolling(1, 2, +1)
case(ichar('2'))
call rolling(1, 2, -1)
case(ichar('x'))
call rolling(0, 1, -1)
case(ichar('y'))
call rolling(0, 1, +1)
case(ichar('w'))
call forwarding(1, +1)
case(ichar('z'))
call forwarding(1, -1)
case(ichar('s'))
call forwarding(0, +1)
case(ichar('a'))
call forwarding(0, -1)
case(ichar(' '))
call forwarding(2, +1)
case(ichar('b'))
call forwarding(2, -1)
case(ichar('1'))
call insordelcube()
case(ichar('m'))
shadowmode= .not. shadowmode
case(ichar('/'))
detailmode= .not. detailmode
case(ichar('r'))
result=messageboxqq('現在の方向を光源に設定しますか？'C, 'Neo cubes'C, MB$OKCANCEL)
if (result .eq. MB$IDOK) then
ray(0:2)=matrix0(2,0:2)
shadowmode=.true.
end if
case(ichar('e'))
result=messageboxqq('視線をリセットしますか？'C, 'Neo cubes'C, MB$OKCANCEL)
if (result .eq. MB$IDOK) call reseteye()
case(ichar('i'))
call loadfile()
case(ichar('o'))
call savefile()
case(ichar('d'))
debugmode=.not.debugmode
case(ichar('h'))
helpmode=.not.helpmode
case(ichar('c'))
screenmode=screenmode+1
if(screenmode .eq. 4) then
screenmode=1
end if
call setcutline()
end select

end do

result=messageboxqq('Qキーが押されました。終了します'C, 'Neo cubes'C, MB$OKCANCEL)
if (result .eq. MB$IDCANCEL) goto 111

end

!メインウインドウの大きさを設定する
logical(4) function initialsettings()
use common_variables
implicit none
integer(4) i
type (qwinfo) qwi
qwi.type=QWIN$MAX !最大化する
i=SetWSizeQQ(QWIN$FRAMEWINDOW, qwi)

initialsettings=.true.
end function initialsettings

!ウインドウの大きさを設定する
subroutine setchildwindowwidth()
use common_variables
implicit none
integer(4) i
type (qwinfo) qwi

qwi.type=QWIN$MIN !最初から表示されているウインドウを最小化する
i=SetWSizeQQ(defaultwindow, qwi)

qwi.type=QWIN$MAX !１を最大化する
i=SetWSizeQQ(childwindow1, qwi)

end subroutine setchildwindowwidth

real(8) function sgn(a)
use common_variables
implicit none

real(8) a

if( a .gt. 0.0_8 ) sgn=1.0_8
if( a .eq. 0.0_8 ) sgn=0.0_8
if( a .lt. 0.0_8 ) sgn=-1.0_8

end function sgn

subroutine setcutline()
use common_variables
implicit none

select case (screenmode)
case(1)
cutline=-0.49_8
cutline2=-0.49_8 !Detail mode用
screen=1.0_8
eyeoffset=1.6_8
zoom=700
case(2)
cutline=-1.9_8
cutline2=-1.9_8 !Detail mode用
screen=1.0_8
eyeoffset=2.0_8 !5.0 または 2.0
zoom=300
case(3)
cutline=0.0_8
cutline2=0.0_8 !Detail mode用
screen=1.0_8
eyeoffset=5.0_8 !5.0 または 2.0
zoom=400
end select

end subroutine setcutline

!matrix0(0:2,0:2)を吸着する
subroutine reseteye()
use common_variables
implicit none

real(8) s, a, s1, a1
integer(2) j, k, k1
real(8), external :: sgn

matrix=0
do j=0, 2
k1=0
s1=matrix0(0,j)
a1=dabs(s1)
s1=sgn(s1)
do k=1, 2
s=matrix0(k,j)
a=dabs(s)
s=sgn(s)
if( a .gt. a1) then
k1=k
a1=a
s1=s
end if
end do
matrix(k1,j)=s1
end do

if( (dabs(matrix(0,0))+dabs(matrix(0,1))+dabs(matrix(0,2)) .eq. 1.0_8) .and. (dabs(matrix(1,0))+dabs(matrix(1,1))+dabs(matrix(1,2)) .eq. 1.0_8) .and. (dabs(matrix(2,0))+dabs(matrix(2,1))+dabs(matrix(2,2)) .eq. 1.0_8) ) then
matrix0=matrix
else
result=messageboxqq('吸着できません。方向(0,0,1)に初期化しますか？'C, 'Neo cubes'C, MB$OKCANCEL)
if (result .eq. MB$IDOK) call reseteye001()
end if

end subroutine reseteye

subroutine reseteye001()
use common_variables
implicit none

matrix0(0,0)=1.0_8
matrix0(0,1)=0.0_8
matrix0(0,2)=0.0_8
matrix0(1,0)=0.0_8
matrix0(1,1)=1.0_8
matrix0(1,2)=0.0_8
matrix0(2,0)=0.0_8
matrix0(2,1)=0.0_8
matrix0(2,2)=1.0_8

end subroutine reseteye001

!視点(eye(0), eye(1), eye(2))を前方に動かす
subroutine forwarding(i, s)
use common_variables
implicit none
integer(4) i, s
real(8) x1(0:2)

x1(0:2)=0.0_8+eye(0:2)

do while( (int4(0.5_8+x1(0)) .eq. eye(0)) .and. (int4(0.5_8+x1(1)) .eq. eye(1)) .and. (int4(0.5_8+x1(2)) .eq. eye(2)))
x1(0:2)=x1(0:2)+s*matrix0(i,0:2)
end do
eye(0:2)=int4(0.5_8+x1(0:2))

end subroutine forwarding

!視点(eye(0),eye(1),eye(2))の前方の立方体を追加／削除する
subroutine insordelcube()
use common_variables
implicit none
real(8) x1(0:2)
integer(4) x2(0:2)
integer(4) j,k

x1(0:2)=0.0_8+eye(0:2)

do while( (int4(0.5_8+x1(0)) .eq. eye(0)) .and. (int4(0.5_8+x1(1)) .eq. eye(1)) .and. (int4(0.5_8+x1(2)) .eq. eye(2)))
x1(0:2)=x1(0:2)+matrix0(2,0:2)
end do
x2(0:2)=int4(0.5_8+x1(0:2))

if( ptrcube .gt. 0 ) then
do j=0, ptrcube-1
if( (x2(0).eq.cube(0,j)) .and. (x2(1).eq.cube(1,j)) .and. (x2(2).eq.cube(2,j)) ) goto 100
end do
end if

if(ptrcube .eq. maxcubes-1) then
result=messageboxqq('これ以上作成できません'C, 'Neo cubes'C, MB$OK)
goto 200
end if

cube(0:2,ptrcube)=x2(0:2)
ptrcube=ptrcube+1
goto 200

100 continue
if(j .lt. ptrcube-1) then
do k=j, ptrcube-2
cube(0:2,k)=cube(0:2,k+1)
end do
end if
ptrcube=ptrcube-1

200 continue

end subroutine insordelcube

!視点の方向 matrix0 を回転する
subroutine rolling(rix, riy, s)
use common_variables
implicit none
integer(2) rix, riy, s
integer(2) j, k
real(8) r

do j=0, 2
r = cos_t * matrix0(rix, j) - s * sin_t * matrix0(riy, j)
matrix0(riy, j) = s * sin_t * matrix0(rix, j) + cos_t * matrix0(riy, j)
matrix0(rix, j) = r
end do

!以下は吸着の処理
do j=0, 2
do k=0, 2
if( (dabs(1-dabs(matrix0(k,j))).gt.0.00001) .and. (dabs(matrix0(k,j)).gt.0.00001) ) goto 999
end do
end do
do j=0, 2
do k=0, 2
if( dabs(1-dabs(matrix0(k,j))).le.0.00001 ) then
if(matrix0(k,j).gt.0) then
r=1.0_8
else
r=-1.0_8
end if
else
r=0.0_8
end if
matrix0(k,j)=r
end do
end do
999 continue

end subroutine rolling

subroutine vram_and_or()
use common_variables
implicit none
vram(rgnxm:rgnxp,rgnym:rgnyp)=( vram(rgnxm:rgnxp,rgnym:rgnyp) .and. vram3(rgnxm:rgnxp,rgnym:rgnyp) )
vram(rgnxm:rgnxp,rgnym:rgnyp)=( vram(rgnxm:rgnxp,rgnym:rgnyp) .or. vram2(rgnxm:rgnxp,rgnym:rgnyp) )
end subroutine vram_and_or

subroutine drawpolygon21(d)
use common_variables
implicit none

integer(2) d
integer(2) x0, y0, x1, y1
integer(2) j4, j5
integer(2) pxx, pyy

plx(0:d-1)=polydata(0:d-1).wx
ply(0:d-1)=polydata(0:d-1).wy
spcbuf=0
dots=d
call poly21()

afterpolydots=spcbuf

if(spcbuf .lt. 3) return !３角形以下では何もしない

!点を整数化し、重複する点を除外する
!すべての点を含む矩形の範囲を知る
j4=0
x0=size_x/2 + int2(cxbuf(ptrcbuf(spcbuf-1)))
rgnxm2=x0
rgnxp2=x0
y0=size_y/2 + int2(cybuf(ptrcbuf(spcbuf-1)))
rgnym2=y0
rgnyp2=y0
do j5=0, spcbuf-1
x1=size_x/2 + int2(cxbuf(ptrcbuf(j5)))
if( x1 .lt. rgnxm2 ) rgnxm2=x1
if( x1 .gt. rgnxp2 ) rgnxp2=x1
y1=size_y/2 + int2(cybuf(ptrcbuf(j5)))
if( y1 .lt. rgnym2 ) rgnym2=y1
if( y1 .gt. rgnyp2 ) rgnyp2=y1
if( (x1 .ne. x0) .or. (y1 .ne. y0) ) then
dispdata(j4).xcoord=x1
dispdata(j4).ycoord=y1
j4=j4+1
x0=x1
y0=y1
end if
end do

afterpolydots=j4

if(j4 .lt. 3) return !３角形以下の場合は何もしない

!リージョンが明らかに重ならないときは、何もしない
if( rgnxp .lt. rgnxm2 ) goto 999
if( rgnyp .lt. rgnym2 ) goto 999
if( rgnxp2 .lt. rgnxm ) goto 999
if( rgnyp2 .lt. rgnym ) goto 999

pxx=dispdata(0).xcoord
pyy=dispdata(0).ycoord

do j5=1, afterpolydots-2
trix(0)=pxx
triy(0)=pyy
trix(1)=dispdata(j5).xcoord
triy(1)=dispdata(j5).ycoord
trix(2)=dispdata(j5+1).xcoord
triy(2)=dispdata(j5+1).ycoord
call drawtriangle21()
end do

return

999 continue
afterpolydots=0

end subroutine drawpolygon21

subroutine drawtriangle21()
use common_variables
implicit none

integer(4) co
integer(2) tmp
real(8) dxybig, dxysmall
real(8) xbig, xsmall
integer(2) left, right
integer(2) j
integer(2) dybig, dysmall
real(8) d

!y座標が小さい順に並べる
if(triy(0) .gt. triy(1)) then
tmp=triy(0)
triy(0)=triy(1)
triy(1)=tmp
tmp=trix(0)
trix(0)=trix(1)
trix(1)=tmp
end if
if(triy(0) .gt. triy(2)) then
tmp=triy(0)
triy(0)=triy(2)
triy(2)=tmp
tmp=trix(0)
trix(0)=trix(2)
trix(2)=tmp
end if
if(triy(1) .gt. triy(2)) then
tmp=triy(2)
triy(2)=triy(1)
triy(1)=tmp
tmp=trix(2)
trix(2)=trix(1)
trix(1)=tmp
end if

dybig=triy(2)-triy(0)
dxybig=(trix(2)-trix(0))/(dybig*2.0_8)

!上半分
dysmall=triy(1)-triy(0)
dxysmall=(trix(1)-trix(0))/(dysmall*2.0_8)
select case (dysmall)

case(0) !横一線
left=imin0( trix(0), trix(1) )
right=imax0( trix(0), trix(1) )
vram2(left:right, triy(0))= ( vram2(left:right, triy(0)) .and. shadowcolor ) !塗る

case default !一般
xbig=trix(0)+dxybig
xsmall=trix(0)+dxysmall
!最初の横一線
left=int2( 0.5_8+ dmin1( trix(0)+0.0_8, xbig, xsmall ) )
right=int2( 0.5_8+ dmax1( trix(0)+0.0_8, xbig, xsmall ) )
vram2(left:right, triy(0))= ( vram2(left:right, triy(0)) .and. shadowcolor ) !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(0)+1, triy(1)-1
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
vram2(left:right, j)= ( vram2(left:right, j) .and. shadowcolor ) !塗る
xbig=xbig+dxybig*2
xsmall=xsmall+dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

end select

!下半分
dysmall=triy(2)-triy(1)
dxysmall=(trix(2)-trix(1))/(dysmall*2.0_8)
select case (dysmall)

case(0) !横一線
left=imin0( trix(1), trix(2) )
right=imax0( trix(1), trix(2) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

case default !一般
xbig=trix(2)-dxybig
xsmall=trix(2)-dxysmall
!最初の横一線
left=int2( 0.5_8+ dmin1( trix(2)+0.0_8, xbig, xsmall ) )
right=int2( 0.5_8+ dmax1( trix(2)+0.0_8, xbig, xsmall ) )
vram2(left:right, triy(2))= ( vram2(left:right, triy(2)) .and. shadowcolor ) !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(2)-1, triy(1)+1, -1
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
vram2(left:right, j)= ( vram2(left:right, j) .and. shadowcolor ) !塗る
xbig=xbig-dxybig*2
xsmall=xsmall-dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

end select

end subroutine drawtriangle21

!(polydata(0).wx,polydata(0).wy) - (polydata(1).wx,polydata(1).wy) - ... -
! (polydata(d-1).wx,polydata(d-1).wy) を色 co で vram にポリゴンを表示する
!afterpolydots=実際に描画された多角形の頂点数
subroutine drawpolygon11(co, d)
use common_variables
implicit none

integer(4) co
integer(2) d
integer(2) x0, y0, x1, y1
integer(2) j4, j5
integer(2) pxx, pyy

plx(0:d-1)=polydata(0:d-1).wx
ply(0:d-1)=polydata(0:d-1).wy
spcbuf=0
dots=d
call poly11()

afterpolydots=spcbuf

if(spcbuf .lt. 3) goto 999 !３角形以下では何もしない

!点を整数化し、重複する点を除外する
!すべての点を含む矩形の範囲を知る
j4=0
x0=size_x/2 + int2(cxbuf(ptrcbuf(spcbuf-1)))
rgnxm=x0
rgnxp=x0
y0=size_y/2 + int2(cybuf(ptrcbuf(spcbuf-1)))
rgnym=y0
rgnyp=y0
do j5=0, spcbuf-1
x1=size_x/2 + int2(cxbuf(ptrcbuf(j5)))
if( x1 .lt. rgnxm ) rgnxm=x1
if( x1 .gt. rgnxp ) rgnxp=x1
y1=size_y/2 + int2(cybuf(ptrcbuf(j5)))
if( y1 .lt. rgnym ) rgnym=y1
if( y1 .gt. rgnyp ) rgnyp=y1
if( (x1 .ne. x0) .or. (y1 .ne. y0) ) then
dispdata(j4).xcoord=x1
dispdata(j4).ycoord=y1
j4=j4+1
x0=x1
y0=y1
end if
end do

afterpolydots=j4

if(j4 .lt. 3) goto 999 !３角形以下の場合は何もしない

vram2(rgnxm:rgnxp, rgnym:rgnyp)=0
vram3(rgnxm:rgnxp, rgnym:rgnyp)=16#ffffff

pxx=dispdata(0).xcoord
pyy=dispdata(0).ycoord

do j5=1, afterpolydots-2
trix(0)=pxx
triy(0)=pyy
trix(1)=dispdata(j5).xcoord
triy(1)=dispdata(j5).ycoord
trix(2)=dispdata(j5+1).xcoord
triy(2)=dispdata(j5+1).ycoord
call drawtriangle11(co)
end do

999 continue

end subroutine drawpolygon11

!整数座標( trix(0), triy(0) ) - ( trix(1), triy(1) ) - ( trix(2), triy(2) )
!の３角形を色 co で vram に表示する
subroutine drawtriangle11(co)
use common_variables
implicit none

integer(4) co
integer(2) tmp
real(8) dxybig, dxysmall
real(8) xbig, xsmall
integer(2) left, right
integer(2) j
integer(2) dybig, dysmall
real(8) d

!y座標が小さい順に並べる
if(triy(0) .gt. triy(1)) then
tmp=triy(0)
triy(0)=triy(1)
triy(1)=tmp
tmp=trix(0)
trix(0)=trix(1)
trix(1)=tmp
end if
if(triy(0) .gt. triy(2)) then
tmp=triy(0)
triy(0)=triy(2)
triy(2)=tmp
tmp=trix(0)
trix(0)=trix(2)
trix(2)=tmp
end if
if(triy(1) .gt. triy(2)) then
tmp=triy(2)
triy(2)=triy(1)
triy(1)=tmp
tmp=trix(2)
trix(2)=trix(1)
trix(1)=tmp
end if

dybig=triy(2)-triy(0)
dxybig=(trix(2)-trix(0))/(dybig*2.0_8)

!上半分
dysmall=triy(1)-triy(0)
dxysmall=(trix(1)-trix(0))/(dysmall*2.0_8)
select case (dysmall)

case(0) !横一線
left=imin0( trix(0), trix(1) )
right=imax0( trix(0), trix(1) )
vram2(left:right, triy(0))= co !塗る
vram3(left:right, triy(0))= 0 !塗る

case default !一般
xbig=trix(0)+dxybig
xsmall=trix(0)+dxysmall
!最初の横一線
left=int2( 0.5_8+ dmin1( trix(0)+0.0_8, xbig, xsmall ) )
right=int2( 0.5_8+ dmax1( trix(0)+0.0_8, xbig, xsmall ) )
vram2(left:right, triy(0))= co !塗る
vram3(left:right, triy(0))= 0 !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(0)+1, triy(1)-1
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
vram2(left:right, j)= co !塗る
vram3(left:right, j)= 0 !塗る
xbig=xbig+dxybig*2
xsmall=xsmall+dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

end select

!下半分
dysmall=triy(2)-triy(1)
dxysmall=(trix(2)-trix(1))/(dysmall*2.0_8)
select case (dysmall)

case(0) !横一線
left=imin0( trix(1), trix(2) )
right=imax0( trix(1), trix(2) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

case default !一般
xbig=trix(2)-dxybig
xsmall=trix(2)-dxysmall
!最初の横一線
left=int2( 0.5_8+ dmin1( trix(2)+0.0_8, xbig, xsmall ) )
right=int2( 0.5_8+ dmax1( trix(2)+0.0_8, xbig, xsmall ) )
vram2(left:right, triy(2))= co !塗る
vram3(left:right, triy(2))= 0 !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(2)-1, triy(1)+1, -1
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
vram2(left:right, j)= co !塗る
vram3(left:right, j)= 0 !塗る
xbig=xbig-dxybig*2
xsmall=xsmall-dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_8+ dmin1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
right=int2( 0.5_8+ dmax1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

end select

end subroutine drawtriangle11

!(plx(0), ply(0)) - (ply(1), ply(1)) - ... - (plx(dots-1), ply(dots-1))
! のポリゴンと四角リージョンの論理積をとる
!cxbuf(spcbuf++), cybuf(spcbuf++) に atan2 でソートされた頂点が入る
!ソートされた結果はポインタの配列 ptrcbuf(0:spcbuf-1) に入っている
subroutine poly11()
use common_variables
implicit none
integer(2), external :: cmp_function_poly11
real(8) poly_a(0:dots-1), poly_b(0:dots-1), poly_c(0:dots-1)
real(8) poly_sign(0:dots-1)
real(8) cx, cy, sign
integer(2) j, k

!四角リージョンの内側にあるポリゴンの頂点は採用する
k=0
do j=0, dots-1
if( plx(j) .ge. 0.0_8-size_x/2 ) then
if( plx(j) .le. size_x/2-1.0_8 ) then
if( ply(j) .ge. 0.0_8-size_y/2 ) then
if( ply(j) .le. size_y/2-1.0_8 ) then
cxbuf(spcbuf)=plx(j)
cybuf(spcbuf)=ply(j)
spcbuf=spcbuf+1
k=k+1
end if
end if
end if
end if
end do

if(k .eq. dots) goto 100 !すべての点が四角リージョンの内側にある場合、すぐにソート処理へ行く

cx=0.0_8
cy=0.0_8
do j=0, dots-1
cx=cx+plx(j)
cy=cy+ply(j)
end do
cx=cx/dots
cy=cy/dots

do j=0, dots-2
line_x0=plx(j)
line_y0=ply(j)
line_x1=plx(j+1)
line_y1=ply(j+1)
call getabc()
poly_a(j)=line_a
poly_b(j)=line_b
poly_c(j)=line_c
end do
line_x0=plx(dots-1)
line_y0=ply(dots-1)
line_x1=plx(0)
line_y1=ply(0)
call getabc()
poly_a(dots-1)=line_a
poly_b(dots-1)=line_b
poly_c(dots-1)=line_c

poly_sign(0:dots-1)= poly_a(0:dots-1)*cx + poly_b(0:dots-1)*cy + poly_c(0:dots-1)

!ポリゴンの内側にある四角リージョンの頂点は採用する
do j=0, 3
do k=0, dots-1
sign= poly_a(k)*rctx(j) + poly_b(k)*rcty(j) + poly_c(k)
if( poly_sign(k)*sign .lt. 0.0_8 ) goto 200 !異符号は採用しない
end do
cxbuf(spcbuf)=rctx(j)
cybuf(spcbuf)=rcty(j)
spcbuf=spcbuf+1
200 continue
end do

!ポリゴンの線と四角リージョンの線の交点が有効なものは採用する
do j=0, dots-2
crx(0)=plx(j)
cry(0)=ply(j)
crx(1)=plx(j+1)
cry(1)=ply(j+1)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()
end do
crx(0)=plx(dots-1)
cry(0)=ply(dots-1)
crx(1)=plx(0)
cry(1)=ply(0)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()

100 continue

!ソートする
if(spcbuf .gt. 0) then
do j=0, spcbuf-1
ptrcbuf(j)=j
end do
end if

if(spcbuf .gt. 3) then

cx=0.0_8
cy=0.0_8
do j=0, spcbuf-1
cx=cx+cxbuf(j)
cy=cy+cybuf(j)
end do
cx=cx/spcbuf
cy=cy/spcbuf
do j=0, spcbuf-1
degree_cbuf(j)= datan2( cybuf(j)-cy, cxbuf(j)-cx )
end do
CALL qsort(ptrcbuf, spcbuf, 2, cmp_function_poly11)

end if

end subroutine poly11

integer(2) function cmp_function_poly11(a1, a2)
use common_variables

implicit none
integer(2) a1, a2
real(8) b1, b2

b1=degree_cbuf(a1)
b2=degree_cbuf(a2)

if( b1 .gt. b2 ) then
cmp_function_poly11=int2(-1)
return
end if

if( b1 .lt. b2 ) then
cmp_function_poly11=int2(1)
return
end if

cmp_function_poly11=int2(0)
return

end function cmp_function_poly11

subroutine poly21()
use common_variables
implicit none
integer(2), external :: cmp_function_poly11
real(8) poly_a(0:dots-1), poly_b(0:dots-1), poly_c(0:dots-1)
real(8) poly_sign(0:dots-1)
real(8) cx, cy, sign
integer(2) j, k
real(8) rgnxmr, rgnxpr, rgnymr, rgnypr, x, y

!四角リージョンの内側にあるポリゴンの頂点は採用する
rgnxmr=plx(0)
rgnxpr=rgnxmr
rgnymr=ply(0)
rgnypr=rgnymr
k=0
do j=0, dots-1
x=plx(j)
if( x .lt. rgnxmr ) rgnxmr=x
if( x .gt. rgnxpr ) rgnxpr=x
y=ply(j)
if( y .lt. rgnymr ) rgnymr=y
if( y .gt. rgnypr ) rgnypr=y

if( x .ge. 0.0_8-size_x/2 ) then
if( x .le. size_x/2-1.0_8 ) then
if( y .ge. 0.0_8-size_y/2 ) then
if( y .le. size_y/2-1.0_8 ) then
cxbuf(spcbuf)=x
cybuf(spcbuf)=y
spcbuf=spcbuf+1
k=k+1
end if
end if
end if
end if
end do

!リージョンが明らかに重ならないときは、何もしない
if( 0.0_8+rgnxp .lt. size_x/2+rgnxmr ) goto 999
if( 0.0_8+rgnyp .lt. size_y/2+rgnymr ) goto 999
if( size_x/2+rgnxpr .lt. 0.0_8+rgnxm ) goto 999
if( size_y/2+rgnypr .lt. 0.0_8+rgnym ) goto 999

if(k .eq. dots) goto 100 !すべての点が四角リージョンの内側にある場合、すぐにソート処理へ行く

cx=0.0_8
cy=0.0_8
do j=0, dots-1
cx=cx+plx(j)
cy=cy+ply(j)
end do
cx=cx/dots
cy=cy/dots

do j=0, dots-2
line_x0=plx(j)
line_y0=ply(j)
line_x1=plx(j+1)
line_y1=ply(j+1)
call getabc()
poly_a(j)=line_a
poly_b(j)=line_b
poly_c(j)=line_c
end do
line_x0=plx(dots-1)
line_y0=ply(dots-1)
line_x1=plx(0)
line_y1=ply(0)
call getabc()
poly_a(dots-1)=line_a
poly_b(dots-1)=line_b
poly_c(dots-1)=line_c

poly_sign(0:dots-1)= poly_a(0:dots-1)*cx + poly_b(0:dots-1)*cy + poly_c(0:dots-1)

!ポリゴンの内側にある四角リージョンの頂点は採用する
do j=0, 3
do k=0, dots-1
sign= poly_a(k)*rctx(j) + poly_b(k)*rcty(j) + poly_c(k)
if( poly_sign(k)*sign .lt. 0.0_8 ) goto 200 !異符号は採用しない
end do
cxbuf(spcbuf)=rctx(j)
cybuf(spcbuf)=rcty(j)
spcbuf=spcbuf+1
200 continue
end do

!ポリゴンの線と四角リージョンの線の交点が有効なものは採用する
do j=0, dots-2
crx(0)=plx(j)
cry(0)=ply(j)
crx(1)=plx(j+1)
cry(1)=ply(j+1)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()
end do
crx(0)=plx(dots-1)
cry(0)=ply(dots-1)
crx(1)=plx(0)
cry(1)=ply(0)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()

100 continue

!ソートする
if(spcbuf .gt. 0) then
do j=0, spcbuf-1
ptrcbuf(j)=j
end do
end if

if(spcbuf .gt. 3) then

cx=0.0_8
cy=0.0_8
do j=0, spcbuf-1
cx=cx+cxbuf(j)
cy=cy+cybuf(j)
end do
cx=cx/spcbuf
cy=cy/spcbuf
do j=0, spcbuf-1
degree_cbuf(j)= datan2( cybuf(j)-cy, cxbuf(j)-cx )
end do
CALL qsort(ptrcbuf, spcbuf, 2, cmp_function_poly11)

end if

return

999 continue
spcbuf=0

end subroutine poly21

subroutine displaydetail()
use common_variables
implicit none

integer(2) j, k, j3, j4, j5, j9
real(8) x0, y0, z0, x1, y1, z1, aa, bb, cc, dd, x2, y2, z2
real(8) xx0(0:7), yy0(0:7), zz0(0:7), t(0:3), dz
integer(2) n, ns
real(8) sign, sign2, t1, sign3
integer(4) j1, k1, j2, j6, j7, k2, j8
integer(4) ix, iy, iz
integer(4) ipiv(3)
integer(2), external :: cmp_function_polygon

!matrix = inverse transpose matrix0 逆行列を求める
matrix=matrix0
call getrf(matrix, ipiv)
call getri(matrix, ipiv)
!matrix を転置する
matrix=transpose(matrix)

!面が何個表示されたか数える（デバッグ用）*******
surfacecount=0

if(ptrcube .eq. 0) goto 999

!表示処理

!cubeoffset(xyz, 頂点番号, 面番号)(定数)→rotoffset(xyz, 頂点番号, 面番号)ベクトル回転する
do j=0, 5
do k=0, 3
rotoffset(0:2,k,j)=matmul(matrix,cubeoffset(0:2,k,j))
end do
end do

do j=0, 5

!rot_abc(面番号) 各面について、a,b,cを計算する
sur_x0=rotoffset(0,0,j)
sur_y0=rotoffset(1,0,j)
sur_z0=rotoffset(2,0,j)
sur_x1=rotoffset(0,1,j)
sur_y1=rotoffset(1,1,j)
sur_z1=rotoffset(2,1,j)
sur_x2=rotoffset(0,2,j)
sur_y2=rotoffset(1,2,j)
sur_z2=rotoffset(2,2,j)
call getabcd()
rot_a(j)=sur_a
rot_b(j)=sur_b
rot_c(j)=sur_c

!signc(面番号) 立方体の真ん中の点について、６つの面に対する正負を計算
!ここでは真ん中の点は(0,0,0)なので
!符号=a*x + b*y + c*z + d よって符号=dになる
!dは直前のgetabcd()でsur_dに入っている
signc(j)= sur_d

end do

!すべての頂点座標を計算する
do j1=0, ptrcube-1
!各立方体の真ん中の点を回転
seev(0:2)=0.0_8+cube(0:2,j1)-eye(0:2)
rot(0:2,j1)=matmul(matrix,seev)

do j=0, 5

!sur(xyz, 頂点番号, 面番号, 通し番号)頂点の座標
do k=0, 3
sur(0:2, k, j, j1) = rot(0:2, j1) + rotoffset(0:2, k, j)
end do

!rot_d(面番号, 通し番号) dを計算する
rot_d(j,j1)= -( rot_a(j)*sur(0,0,j,j1) + rot_b(j)*sur(1,0,j,j1) + rot_c(j)*sur(2,0,j,j1) )

!actcontact(面番号, 通し番号)互いに接している面は除外する
ix=cube(0,j1)+dircontoffset(0,j)
iy=cube(1,j1)+dircontoffset(1,j)
iz=cube(2,j1)+dircontoffset(2,j)
do k1=0, ptrcube-1
if( cube(0,k1).eq.ix ) then
if( cube(1,k1).eq.iy ) then
if( cube(2,k1).eq.iz ) then
actcontact(j,j1)=.false. !除外する
goto 100
end if
end if
end if
end do
actcontact(j,j1)=.true. !採用する
100 continue

end do
end do

endofsortptr=0
do j1=0, ptrcube-1

!sur(xyz, 頂点番号, 面番号, 通し番号)
!１ヶの立方体の座標の中で、どれかがカットラインより前に有れば採用する可能性がある
do j=0, 5
do k=0, 3
if( sur(2,k,j,j1) .gt. cutline2 ) goto 210
end do
end do
!すべての座標がカットラインより後ろにあった
do j=0, 5
act(j,j1)=.false.
end do
goto 200

210 continue

!ソート用ポインタに登録
sortptr(endofsortptr)=j1
endofsortptr=endofsortptr+1

!act(面番号, 通し番号)
!actcontactがfalseのときは無条件でfalse
!視点(0,0,-eyeoffset)に対して隠れている３つの面を除外する
do j=0, 5
act(j,j1)=actcontact(j,j1)
if( signc(j) * ( rot_c(j)*-eyeoffset + rot_d(j,j1) ) .gt. 0.0_8 ) act(j,j1)=.false. !正負が合っていれば、隠れ面
end do

200 continue
end do

!影が何個表示されたか数える（デバッグ用）*******
shadowcount=0

!影モードの準備処理（表示する立方体が無いときは、飛ばして良し）
if( endofsortptr .gt. 0 ) then

!光源を回転
rotray(0:2)=matmul(matrix,ray(0:2))

do j=0, 5
!actshadow(面番号)
actshadow(j)=.true.
!rotoffset(xyz, 頂点番号, 面番号)
!各面の真ん中の点から光源の0.1倍のベクトルを伸ばした座標を計算する
x0=( rotoffset(0,0,j)+rotoffset(0,1,j)+rotoffset(0,2,j)+rotoffset(0,3,j) ) / 4 + rotray(0)/10
y0=( rotoffset(1,0,j)+rotoffset(1,1,j)+rotoffset(1,2,j)+rotoffset(1,3,j) ) / 4 + rotray(1)/10
z0=( rotoffset(2,0,j)+rotoffset(2,1,j)+rotoffset(2,2,j)+rotoffset(2,3,j) ) / 4 + rotray(2)/10
!この立方体の真ん中の点(0,0,0)に対する各面の正／負と比較し、
!合っていなければtrueの可能性有り、合っていればfalse
do k=0, 5
sign= rot_a(k)*x0 + rot_b(k)*y0 + rot_c(k)*z0 + signc(k) !"rot_d(k)"はsignc(k)と同じである
if( (dabs(sign) .gt. 0.00001_8) .and. (sign * signc(k) .lt. 0.0_8) ) goto 300 !合っていないので採用される可能性ある
end do
!全部合っているということは、光源ベクトルの方向は、この面から立方体の内側へ向かっている
actshadow(j)=.false.
300 continue
end do

!centerpoint(xyz, 面番号, 通し番号) 各面の真ん中の点
!sur(xyz, 頂点番号, 面番号, 通し番号)の各頂点を足して４で割る
centerpoint(0:2, 0:5, 0:ptrcube-1)= ( sur(0:2, 0, 0:5, 0:ptrcube-1)+sur(0:2, 1, 0:5, 0:ptrcube-1)+sur(0:2, 2, 0:5, 0:ptrcube-1)+sur(0:2, 3, 0:5, 0:ptrcube-1) )/4

end if

!ポリゴンを表示する
if( endofsortptr .gt. 0 ) then

!視点(0,0,-eyeoffset)からの距離を計算
distance(sortptr(0:endofsortptr-1))= rot(0,sortptr(0:endofsortptr-1))**2 + rot(1,sortptr(0:endofsortptr-1))**2 + (rot(2,sortptr(0:endofsortptr-1))+eyeoffset)**2
!ソートする
CALL qsort(sortptr, endofsortptr, 4, cmp_function_polygon)

!表示する
do j1=0, endofsortptr-1
k1=sortptr(j1)
do j=0, 5
if( act(j,k1) ) then

xx0(0:3)=sur(0, 0:3, j, k1)
yy0(0:3)=sur(1, 0:3, j, k1)
zz0(0:3)=sur(2, 0:3, j, k1)
if( (zz0(0) .le. cutline2) .and. (zz0(1) .le. cutline2) .and. (zz0(2) .le. cutline2) .and. (zz0(3) .le. cutline2) ) then
goto 310 !すべての頂点がカットラインの後ろにあるため、何もしない
end if
!どれかの頂点がカットラインより前にある
xx0(4:7)=xx0(0:3)
yy0(4:7)=yy0(0:3)
zz0(4:7)=zz0(0:3)
do j4=0, 3
if(zz0(j4) .le. cutline2 ) then
!カットラインより後ろにある
!どれかは前にあり、どれかは後ろにあるので、３次元座標を切り取る
if(zz0(0) .le. cutline2 ) then
j9=0
do while( zz0(j9) .le. cutline2 )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline2 ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
else
j5=0
do while( zz0(j5) .gt. cutline2 )
j5=j5+1
end do
j9=j5
do while( zz0(j9) .le. cutline2 )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline2 ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
end if
!zz0(j9-1)とzz0(j9)の切り取り点, 点zz0(j9)～zz0(j5-1), zz0(j5-1)とzz0(j5)の切り取り点

dz= zz0(j9)-zz0(j9-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 310 !zz0(j9-1)とzz0(j9)が近い場合、何もしない
dz= ( cutline2-zz0(j9-1) ) / dz
xx0(j9-1)= xx0(j9-1) + ( xx0(j9)-xx0(j9-1) ) * dz
yy0(j9-1)= yy0(j9-1) + ( yy0(j9)-yy0(j9-1) ) * dz
zz0(j9-1)= 0.0_8+cutline2

dz= zz0(j5)-zz0(j5-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 310 !zz0(j5-1)とzz0(j5)が近い場合、何もしない
dz= ( cutline2-zz0(j5-1) ) / dz
xx0(j5)= xx0(j5-1) + ( xx0(j5)-xx0(j5-1) ) * dz
yy0(j5)= yy0(j5-1) + ( yy0(j5)-yy0(j5-1) ) * dz
zz0(j5)= 0.0_8+cutline2

!投影する
zz0(j9-1:j5) = (screen+eyeoffset)/(zz0(j9-1:j5)+eyeoffset)
ns= j5 - (j9-1) + 1
polydata(0:ns-1).wx = xx0(j9-1:j5)*zz0(j9-1:j5)*zoom
polydata(0:ns-1).wy = yy0(j9-1:j5)*zz0(j9-1:j5)*zoom
goto 320
end if
end do
!４点すべてカットラインより前にあった
!投影する
zz0(0:3) = (screen+eyeoffset)/(zz0(0:3)+eyeoffset)
polydata(0:3).wx = xx0(0:3)*zz0(0:3)*zoom
polydata(0:3).wy = yy0(0:3)*zz0(0:3)*zoom
ns=4

320 continue

surfacecount=surfacecount+1

!影の表示
aa= rot_a(j)
bb= rot_b(j)
cc= rot_c(j)
dd= rot_d(j, k1)
!この面番号jの面に対し、視点(0,0,-eyeoffset)の正／負を調べる
sign= cc*-eyeoffset + dd
!このjの面の頂点番号0から、光源ベクトルを伸ばした座標を計算する
x0= sur(0, 0, j, k1) + rotray(0)
y0= sur(1, 0, j, k1) + rotray(1)
z0= sur(2, 0, j, k1) + rotray(2)
!その正／負を調べる
sign2= aa*x0 + bb*y0 + cc*z0 + dd
!これが０に近い（面と光源が平行）とき、この面は全部影にする
!また、符号が合っている場合、光源に対して面の裏側を見ていることになるので、この面は全部影になる
if( (dabs(sign2) .lt. 0.00001_8) .or. (sign * sign2 .ge. 0.0_8) ) then
!面と光源が平行であるか、符号が合っている
call drawpolygon11( cubecolor(j) .and. shadowcolor, ns ) !面番号に対応した影の色にする
if( afterpolydots .ge. 3 ) then
call vram_and_or()
shadowcount=shadowcount+1
end if
goto 400 !この面の処理終わり
else
!合っていない場合、以下の処理によってこの面は影で上書きされる
call drawpolygon11( cubecolor(j), ns ) !面番号に対応した影でない色にする
if(afterpolydots .lt. 3) goto 400 !影で上書きされる予定のjの面が実際に表示されなかった場合、この面の処理終わり
end if

!分母が0に近いとき、影の頂点は計算不能
t1 = -aa*rotray(0)-bb*rotray(1)-cc*rotray(2)
if( dabs(t1) .lt. 0.00001 ) goto 410

!登録されたすべての立方体の６面が影の発生原因になる可能性ある
do j2=0, ptrcube-1
!通し番号が自分自身でないものについて処理する
if(j2.ne.k1) then
!このj2の立方体の真ん中の点に対する符号を見る
!rot(xyz, 通し番号) 回転後の真ん中の点
sign= aa*rot(0,j2) + bb*rot(1,j2) + cc*rot(2,j2) + dd
!このjの面の頂点番号0から、光源ベクトルを伸ばした座標の正／負はsign2に入っている
if(sign * sign2 .ge. 0.0_8 ) goto 500 !符号が合っているときは影の原因にならないので除外する
do j3=0, 5
!actshadow(面番号)がtrueのものについて処理する
!actcontact(面番号, 通し番号)がfalseのときは無条件でfalse
if( actshadow(j3) .and. actcontact(j3,j2) ) then
!このj3の面の真ん中の点の符号を見る
!centerpoint(xyz, 面番号, 通し番号)
sign= aa*centerpoint(0,j3,j2) + bb*centerpoint(1,j3,j2) + cc*centerpoint(2,j3,j2) + dd
if( dabs(sign) .lt. 0.00001 ) goto 510 !jの面の延長上にあるこの面は影の原因にならない

!j,k1の面の頂点とj2,j3の面の頂点が一致しているか調べる
j4=0
do j5=0, 3
do j9=0, 3
if( dabs(sur(0, j5, j, k1) - sur(0, j9, j3, j2)) + dabs(sur(1, j5, j, k1) - sur(1, j9, j3, j2)) + dabs(sur(2, j5, j, k1) - sur(2, j9, j3, j2)) .lt. 0.00001 ) then
j4=j4+1
x2=sur(0,j5,j,k1)
y2=sur(1,j5,j,k1)
z2=sur(2,j5,j,k1)
end if
end do
end do
if(j4 .ne. 1) goto 600 !1点で接している場合以外は、見逃す

!1点で接している場合、「黒点」の問題が存在する可能性がある
!接している点は(x2, y2, z2)に入っている
!この点から光源ベクトル/10を伸ばす
x2=x2+rotray(0)/10
y2=y2+rotray(1)/10
z2=z2+rotray(2)/10
do j4=0, 5
sign3= rot_a(j4)*x2 + rot_b(j4)*y2 + rot_c(j4)*z2 + rot_d(j4,k1)
if( sign3*signc(j4) .lt. 0.0_8 ) goto 510 !k1の立方体の外側へ向かっているので、これは影の原因ではない
end do

600 continue
!次の方程式をtについて解く
!x == x0 + (x1 - x0)*t
!y == y0 + (y1 - y0)*t
!z == z0 + (z1 - z0)*t
!a*x + b*y + c*z + d == 0
!結果は
!t = (d + a x0 + b y0 + c z0)
! / (a x0 - a x1 + b y0 - b y1 + c z0 - c z1)
!sur(xyz, 頂点番号, 面番号, 通し番号)
xx0(0:3)=sur(0, 0:3, j3, j2)
yy0(0:3)=sur(1, 0:3, j3, j2)
zz0(0:3)=sur(2, 0:3, j3, j2)
t(0:3) = (dd + aa*xx0(0:3) + bb*yy0(0:3) + cc*zz0(0:3)) / t1
!投影する
zz0(0:3)=zz0(0:3)+rotray(2)*t(0:3)
do j4=0, 3
if(zz0(j4) .le. cutline2 ) then
!カットラインより後ろにある
if( dmax1(zz0(0), zz0(1), zz0(2), zz0(3)) .le. cutline2 ) goto 510 !全部後ろにある場合、何もしない
!どれかは前にあり、どれかは後ろにあるので、影の３次元座標を切り取る
zz0(4:7)=zz0(0:3)
xx0(0:3)=xx0(0:3)+rotray(0)*t(0:3)
yy0(0:3)=yy0(0:3)+rotray(1)*t(0:3)
xx0(4:7)=xx0(0:3)
yy0(4:7)=yy0(0:3)
if(zz0(0) .le. cutline2 ) then
j9=0
do while( zz0(j9) .le. cutline2 )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline2 ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
else
j5=0
do while( zz0(j5) .gt. cutline2 )
j5=j5+1
end do
j9=j5
do while( zz0(j9) .le. cutline2 )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline2 ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
end if
!zz0(j9-1)とzz0(j9)の切り取り点, 点zz0(j9)～zz0(j5-1), zz0(j5-1)とzz0(j5)の切り取り点

dz= zz0(j9)-zz0(j9-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 510 !zz0(j9-1)とzz0(j9)が近い場合、何もしない
dz= ( cutline2-zz0(j9-1) ) / dz
xx0(j9-1)= xx0(j9-1) + ( xx0(j9)-xx0(j9-1) ) * dz
yy0(j9-1)= yy0(j9-1) + ( yy0(j9)-yy0(j9-1) ) * dz
zz0(j9-1)= 0.0_8+cutline2

dz= zz0(j5)-zz0(j5-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 510 !zz0(j5-1)とzz0(j5)が近い場合、何もしない
dz= ( cutline2-zz0(j5-1) ) / dz
xx0(j5)= xx0(j5-1) + ( xx0(j5)-xx0(j5-1) ) * dz
yy0(j5)= yy0(j5-1) + ( yy0(j5)-yy0(j5-1) ) * dz
zz0(j5)= 0.0_8+cutline2

zz0(j9-1:j5) = (screen+eyeoffset)/(zz0(j9-1:j5)+eyeoffset)
n= j5 - (j9-1) + 1
polydata(0:n-1).wx = xx0(j9-1:j5)*zz0(j9-1:j5)*zoom
polydata(0:n-1).wy = yy0(j9-1:j5)*zz0(j9-1:j5)*zoom
goto 430
end if
end do
!４点すべてカットラインより前にあった
zz0(0:3) = (screen+eyeoffset)/(zz0(0:3)+eyeoffset)
polydata(0:3).wx = (xx0(0:3) + rotray(0)*t(0:3))*zz0(0:3)*zoom
polydata(0:3).wy = (yy0(0:3) + rotray(1)*t(0:3))*zz0(0:3)*zoom
n=4

430 continue

call drawpolygon21( n )
if( afterpolydots .ge. 3 ) shadowcount=shadowcount+1 !３角形以上のとき、採用された
end if
510 continue
end do
500 continue
end if
end do
!影の原因になる可能性のあるすべての立方体を調べた
410 continue
call vram_and_or()

400 continue

end if
310 continue

end do
end do

end if

999 continue

end subroutine displaydetail

subroutine displaycube()
use common_variables
implicit none

integer(2) j, k, j3, j4, j5, j9
real(8) x0, y0, z0, x1, y1, z1, aa, bb, cc, dd
real(8) xx0(0:7), yy0(0:7), zz0(0:7), t(0:3), dz
integer(2) n, ns
real(8) sign, sign2, t1
integer(4) j1, k1, j2, j6, j7, k2, j8, l
integer(4) oldk1
integer(4) ptrcapture, beforeptrcapture, nextptrcapture, dx, dy, c, p, p1
integer(4) ix, iy, iz
integer(4) ipiv(3)
logical(2) reg
integer(2), external :: cmp_function_polygon

!matrix = inverse transpose matrix0 逆行列を求める
matrix=matrix0
call getrf(matrix, ipiv)
call getri(matrix, ipiv)
!matrix を転置する
matrix=transpose(matrix)

!面が何個表示されたか数える（デバッグ用）*******
surfacecount=0

ptrcapture=0
endofcutangles=0

if(ptrcube .eq. 0) goto 999

!表示処理

!cubeoffset(xyz, 頂点番号, 面番号)(定数)→rotoffset(xyz, 頂点番号, 面番号)ベクトル回転する
do j=0, 5
do k=0, 3
rotoffset(0:2,k,j)=matmul(matrix,cubeoffset(0:2,k,j))
end do
end do

do j=0, 5

!rot_abc(面番号) 各面について、a,b,cを計算する
sur_x0=rotoffset(0,0,j)
sur_y0=rotoffset(1,0,j)
sur_z0=rotoffset(2,0,j)
sur_x1=rotoffset(0,1,j)
sur_y1=rotoffset(1,1,j)
sur_z1=rotoffset(2,1,j)
sur_x2=rotoffset(0,2,j)
sur_y2=rotoffset(1,2,j)
sur_z2=rotoffset(2,2,j)
call getabcd()
rot_a(j)=sur_a
rot_b(j)=sur_b
rot_c(j)=sur_c

!signc(面番号) 立方体の真ん中の点について、６つの面に対する正負を計算
!ここでは真ん中の点は(0,0,0)なので
!符号=a*x + b*y + c*z + d よって符号=dになる
!dは直前のgetabcd()でsur_dに入っている
signc(j)= sur_d

end do

!すべての頂点座標を計算する
do j1=0, ptrcube-1
!各立方体の真ん中の点を回転
seev(0:2)=0.0_8+cube(0:2,j1)-eye(0:2)
rot(0:2,j1)=matmul(matrix,seev)

do j=0, 5

!sur(xyz, 頂点番号, 面番号, 通し番号)頂点の座標
do k=0, 3
sur(0:2, k, j, j1) = rot(0:2, j1) + rotoffset(0:2, k, j)
end do

!rot_d(面番号, 通し番号) dを計算する
rot_d(j,j1)= -( rot_a(j)*sur(0,0,j,j1) + rot_b(j)*sur(1,0,j,j1) + rot_c(j)*sur(2,0,j,j1) )

!actcontact(面番号, 通し番号)互いに接している面は除外する
ix=cube(0,j1)+dircontoffset(0,j)
iy=cube(1,j1)+dircontoffset(1,j)
iz=cube(2,j1)+dircontoffset(2,j)
do k1=0, ptrcube-1
if( cube(0,k1).eq.ix ) then
if( cube(1,k1).eq.iy ) then
if( cube(2,k1).eq.iz ) then
actcontact(j,j1)=.false. !除外する
goto 100
end if
end if
end if
end do
actcontact(j,j1)=.true. !採用する
100 continue

end do
end do

endofcutangles=0
endofsortptr=0
do j1=0, ptrcube-1

if( rot(2,j1) .gt. cutline+sqr3d2 ) then
!立方体の中心が cutline+sqr(3)/2 よりも前にある場合、登録する
!act(面番号, 通し番号)
!actcontactがfalseのときは無条件でfalse
!視点(0,0,-eyeoffset)に対して隠れている３つの面を除外する
do j=0, 5
act(j,j1)=actcontact(j,j1)
if( signc(j) * ( rot_c(j)*-eyeoffset + rot_d(j,j1) ) .gt. 0.0_8 ) act(j,j1)=.false. !正負が合っていれば、隠れ面
end do
ptrcut(0:5, j1)=-1 !カットライン処理とは無関係のため、ポインタをリセットする
goto 201
end if

if( rot(2,j1) .le. cutline-sqr3d2 ) goto 200 !立方体の中心が cutline-sqr(3)/2 よりも後ろにある場合、登録しない

!立方体の中心がカットライン付近にある場合
!act(面番号, 通し番号)
!actcontactがfalseのときは無条件でfalse
!視点(0,0,-eyeoffset)に対して隠れている３つの面を除外する
do j=0, 5
act(j,j1)=actcontact(j,j1)
if( signc(j) * ( rot_c(j)*-eyeoffset + rot_d(j,j1) ) .gt. 0.0_8 ) act(j,j1)=.false. !正負が合っていれば、隠れ面
end do

reg=.false.
do j=0, 5
if( act(j,j1) ) then
xx0(0:3)=sur(0, 0:3, j, j1)
yy0(0:3)=sur(1, 0:3, j, j1)
zz0(0:3)=sur(2, 0:3, j, j1)
if( (zz0(0) .le. cutline) .and. (zz0(1) .le. cutline) .and. (zz0(2) .le. cutline) .and. (zz0(3) .le. cutline) ) then
act(j,j1)= .false. !すべての頂点がカットラインの後ろにあるため、登録しない
goto 210
end if
if( (zz0(0) .gt. cutline) .and. (zz0(1) .gt. cutline) .and. (zz0(2) .gt. cutline) .and. (zz0(3) .gt. cutline) ) then
!すべての頂点がカットラインの前にあるため、登録する
!ポインタをリセットし、通常の４角形を投影させる
ptrcut(j, j1)=-1
reg=.true.
goto 210
end if
!どれかの頂点がカットラインより前にあり、どれかは後ろにある
xx0(4:7)=xx0(0:3)
yy0(4:7)=yy0(0:3)
zz0(4:7)=zz0(0:3)
do j4=0, 3
if(zz0(j4) .le. cutline ) then
!カットラインより後ろにある
!どれかは前にあり、どれかは後ろにあるので、３次元座標を切り取る
if(zz0(0) .le. cutline ) then
j9=0
do while( zz0(j9) .le. cutline )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
else
j5=0
do while( zz0(j5) .gt. cutline )
j5=j5+1
end do
j9=j5
do while( zz0(j9) .le. cutline )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
end if
!zz0(j9-1)とzz0(j9)の切り取り点, 点zz0(j9)～zz0(j5-1), zz0(j5-1)とzz0(j5)の切り取り点

dz= zz0(j9)-zz0(j9-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 200 !zz0(j9-1)とzz0(j9)が近い場合、何もしない
dz= ( cutline-zz0(j9-1) ) / dz
xx0(j9-1)= xx0(j9-1) + ( xx0(j9)-xx0(j9-1) ) * dz
yy0(j9-1)= yy0(j9-1) + ( yy0(j9)-yy0(j9-1) ) * dz
zz0(j9-1)= 0.0_8+cutline

dz= zz0(j5)-zz0(j5-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 200 !zz0(j5-1)とzz0(j5)が近い場合、何もしない
dz= ( cutline-zz0(j5-1) ) / dz
xx0(j5)= xx0(j5-1) + ( xx0(j5)-xx0(j5-1) ) * dz
yy0(j5)= yy0(j5-1) + ( yy0(j5)-yy0(j5-1) ) * dz
zz0(j5)= 0.0_8+cutline

!投影する
zz0(j9-1:j5) = (screen+eyeoffset)/(zz0(j9-1:j5)+eyeoffset)
ns= j5 - (j9-1) + 1
polydata(0:ns-1).wx = xx0(j9-1:j5)*zz0(j9-1:j5)*zoom
polydata(0:ns-1).wy = yy0(j9-1:j5)*zz0(j9-1:j5)*zoom

if(endofcutangles .ge. maxcutangles-ns*2-1) then
result=messageboxqq('カットライン処理領域が足りません'C, 'Neo cubes'C, MB$OK)
goto 999
end if
!ポインタを登録
ptrcut(j, j1)=endofcutangles
cutangles(endofcutangles)=ns+0.0_8
do k=0, ns-1
cutangles(endofcutangles+1+k*2+0)=polydata(k).wx
cutangles(endofcutangles+1+k*2+1)=polydata(k).wy
end do
endofcutangles=endofcutangles+1+ns*2
reg=.true.
goto 210
end if
end do
!ここはバグ
result=messageboxqq('BUG カットライン処理'C, 'Neo cubes'C, MB$OK)

210 continue
end if
end do
if( .not. reg ) goto 200 !0～5の面番号すべてで登録されない場合、ソートポインタにも登録しない

201 continue
!ソート用ポインタに登録
sortptr(endofsortptr)=j1
endofsortptr=endofsortptr+1

200 continue
end do

!影が何個表示されたか数える（デバッグ用）*******
shadowcount=0

!影モードの準備処理（表示する立方体が無いときは、飛ばして良し）
if( (endofsortptr .gt. 0) .and. shadowmode) then

!光源を回転
rotray(0:2)=matmul(matrix,ray(0:2))

do j=0, 5
!actshadow(面番号)
actshadow(j)=.true.
!rotoffset(xyz, 頂点番号, 面番号)
!各面の真ん中の点から光源の0.1倍のベクトルを伸ばした座標を計算する
x0=( rotoffset(0,0,j)+rotoffset(0,1,j)+rotoffset(0,2,j)+rotoffset(0,3,j) ) / 4 + rotray(0)/10
y0=( rotoffset(1,0,j)+rotoffset(1,1,j)+rotoffset(1,2,j)+rotoffset(1,3,j) ) / 4 + rotray(1)/10
z0=( rotoffset(2,0,j)+rotoffset(2,1,j)+rotoffset(2,2,j)+rotoffset(2,3,j) ) / 4 + rotray(2)/10
!この立方体の真ん中の点(0,0,0)に対する各面の正／負と比較し、
!合っていなければtrueの可能性有り、合っていればfalse
do k=0, 5
sign= rot_a(k)*x0 + rot_b(k)*y0 + rot_c(k)*z0 + signc(k) !"rot_d(k)"はsignc(k)と同じである
if( (dabs(sign) .gt. 0.00001_8) .and. (sign * signc(k) .lt. 0.0_8) ) goto 300 !合っていないので採用される可能性ある
end do
!全部合っているということは、光源ベクトルの方向は、この面から立方体の内側へ向かっている
actshadow(j)=.false.
300 continue
end do

!centerpoint(xyz, 面番号, 通し番号) 各面の真ん中の点
!sur(xyz, 頂点番号, 面番号, 通し番号)の各頂点を足して４で割る
centerpoint(0:2, 0:5, 0:ptrcube-1)= ( sur(0:2, 0, 0:5, 0:ptrcube-1)+sur(0:2, 1, 0:5, 0:ptrcube-1)+sur(0:2, 2, 0:5, 0:ptrcube-1)+sur(0:2, 3, 0:5, 0:ptrcube-1) )/4

end if

!ポリゴンを表示する
if( endofsortptr .gt. 0 ) then

!視点(0,0,-eyeoffset)からの距離を計算
distance(sortptr(0:endofsortptr-1))= rot(0,sortptr(0:endofsortptr-1))**2 + rot(1,sortptr(0:endofsortptr-1))**2 + (rot(2,sortptr(0:endofsortptr-1))+eyeoffset)**2
!ソートする
CALL qsort(sortptr, endofsortptr, 4, cmp_function_polygon)

!同じdを持つ平面をグループにする
ppgroup=0
pgroupforcube=0
ppsort=0 !キャプチャー領域を指す予定のsortc(maxcubes*6)へのポインタ（ソートされた順に並んでいる）
do j1=0, endofsortptr-1
k1=sortptr(j1)
do j=0, 5
if( act(j,k1) ) then
dd= rot_d(j,k1)
j2=0
do while (j2 .lt. ppgroup)
if( j .eq. groupforsurface(j2) ) then
if( dabs(dd-groupd(j2)) .lt. 0.00001_8 ) then
!同じ面番号で同じdを持つ
j6=endofnext(j2)
endofnext(j2)=pgroupforcube
nextptr(j6)=pgroupforcube
nextptr(pgroupforcube)=-1
groupforserial(pgroupforcube)=k1 !通し番号
sortc(ppsort)=pgroupforcube
ppsort=ppsort+1
pgroupforcube=pgroupforcube+1
goto 350
end if
end if
j2=j2+1
end do
!新しいdの定義
groupd(ppgroup)=dd
ptrgroup(ppgroup)=pgroupforcube
groupforserial(pgroupforcube)=k1 !通し番号
endofnext(ppgroup)=pgroupforcube
nextptr(pgroupforcube)=-1
sortc(ppsort)=pgroupforcube
ppsort=ppsort+1
pgroupforcube=pgroupforcube+1
groupforsurface(ppgroup)=j !面番号
ppgroup=ppgroup+1
end if
350 continue
end do
end do

!表示する
if( ppgroup .gt. 0 ) then
ptrcapture=0
do j1=0, ppgroup-1

identity(0:ptrcube-1)=.true. !このグループに対する影の処理を制限するための変数

j=groupforsurface(j1) !面番号
beforeptrcapture=ptrcapture
c=0
isrgnsetted=.false.
p=ptrgroup(j1)
j2=0
do while (p .ne. -1)
k1=groupforserial(p) !通し番号
if(j2 .eq. 0) oldk1=k1 !このグループ全体を代表して最初の面の通し番号を取る

identity(k1)=.false. !このグループに所属するすべての面は影の発生原因としない

surfacecount=surfacecount+1

p1=ptrcut(j,k1)
if( p1 .eq. -1 ) then
!元の４角形を採用する
!投影する
!sur(xyz, 頂点番号, 面番号, 通し番号)
do k=0, 3
z0=(screen+eyeoffset)/(sur(2, k, j, k1)+eyeoffset)
polydata(k).wx=sur(0, k, j, k1)*z0*zoom
polydata(k).wy=sur(1, k, j, k1)*z0*zoom
end do
call drawpolygon( cubecolor(j), 4 ) !面番号に対応した影でない色にする
else
!カットライン付近で切り取られた面を採用する
ns=int2(cutangles(p1))
do k=0, ns-1
polydata(k).wx=cutangles(p1+1+k*2+0)
polydata(k).wy=cutangles(p1+1+k*2+1)
end do
call drawpolygon( cubecolor(j), ns ) !面番号に対応した影でない色にする
end if

if( afterpolydots .ge. 3 ) then
c=1 !３角形以上が表示された
dx=rgnxp-rgnxm+1
dy=rgnyp-rgnym+1
nextptrcapture= ptrcapture + 4 + dx*dy*2 !キャプチャー領域のサイズはリージョンのサイズ*2とする
if( nextptrcapture .ge. maxcapture ) then
result=messageboxqq('キャプチャー領域が足りません'C, 'Neo cubes'C, MB$OK)
goto 999
end if
groupforserial(p)=ptrcapture !通し番号へのポインタをキャプチャー領域へのポインタに書き換える
capture(ptrcapture+0)=rgnxm
capture(ptrcapture+1)=rgnym
capture(ptrcapture+2)=dx
capture(ptrcapture+3)=dy
ptrcapture=ptrcapture+4
do j3=rgnym, rgnyp
capture(ptrcapture:ptrcapture+dx-1)=vram3(rgnxm:rgnxm+dx-1,j3) !ここではビットマスクを転送する
ptrcapture=ptrcapture+dx
end do
ptrcapture=nextptrcapture
else
groupforserial(p)=-1 !何も表示しないときは、ポインタを無効にする
end if
p=nextptr(p)
j2=1
end do

if( shadowmode .and. (c .eq. 1) ) then
!３角形が１度以上表示されたことがある場合、
!影を表示する

k1=oldk1
aa= rot_a(j)
bb= rot_b(j)
cc= rot_c(j)
dd= rot_d(j, k1)
!このグループの面に対し、視点(0,0,-eyeoffset)の正／負を調べる
sign= cc*-eyeoffset + dd
!この面の頂点番号0から、光源ベクトルを伸ばした座標を計算する
x0= sur(0, 0, j, k1) + rotray(0)
y0= sur(1, 0, j, k1) + rotray(1)
z0= sur(2, 0, j, k1) + rotray(2)
!その正／負を調べる
sign2= aa*x0 + bb*y0 + cc*z0 + dd
!これが０に近い（面と光源が平行）とき、この面は全部影にする
!また、符号が合っている場合、光源に対して面の裏側を見ていることになるので、この面は全部影になる
if( (dabs(sign2) .lt. 0.00001_8) .or. (sign * sign2 .ge. 0.0_8) ) then
!面と光源が平行であるか、符号が合っている
if( isrgnsetted ) then
vram2(rgnxmgroup:rgnxpgroup,rgnymgroup:rgnypgroup)=( vram2(rgnxmgroup:rgnxpgroup,rgnymgroup:rgnypgroup) .and. shadowcolor )
shadowcount=shadowcount+1
end if
goto 400 !この面の処理終わり
end if

!合っていない場合、以下の処理によってこの面は影で上書きされる

!分母が0に近いとき、影の頂点は計算不能
t1 = -aa*rotray(0)-bb*rotray(1)-cc*rotray(2)
if( dabs(t1) .lt. 0.00001 ) goto 400

!登録されたすべての立方体の６面が影の発生原因になる可能性ある
do j2=0, ptrcube-1
!通し番号がこのグループ自身でないものについて処理する
if( identity(j2) ) then
!このj2の立方体の真ん中の点に対する符号を見る
!rot(xyz, 通し番号) 回転後の真ん中の点
sign= aa*rot(0,j2) + bb*rot(1,j2) + cc*rot(2,j2) + dd
!このjの面の頂点番号0から、光源ベクトルを伸ばした座標の正／負はsign2に入っている
if(sign * sign2 .ge. 0.0_8 ) goto 500 !符号が合っているときは影の原因にならないので除外する
do j3=0, 5
!actshadow(面番号)がtrueのものについて処理する
!actcontact(面番号, 通し番号)がfalseのときは無条件でfalse
if( actshadow(j3) .and. actcontact(j3,j2) ) then
!このj3の面の真ん中の点の符号を見る
!centerpoint(xyz, 面番号, 通し番号)
sign= aa*centerpoint(0,j3,j2) + bb*centerpoint(1,j3,j2) + cc*centerpoint(2,j3,j2) + dd
if( dabs(sign) .lt. 0.00001 ) goto 510 !jの面の延長上にあるこの面は影の原因にならない

!次の方程式をtについて解く
!x == x0 + (x1 - x0)*t
!y == y0 + (y1 - y0)*t
!z == z0 + (z1 - z0)*t
!a*x + b*y + c*z + d == 0
!結果は
!t = (d + a x0 + b y0 + c z0)
! / (a x0 - a x1 + b y0 - b y1 + c z0 - c z1)
!sur(xyz, 頂点番号, 面番号, 通し番号)
xx0(0:3)=sur(0, 0:3, j3, j2)
yy0(0:3)=sur(1, 0:3, j3, j2)
zz0(0:3)=sur(2, 0:3, j3, j2)
t(0:3) = (dd + aa*xx0(0:3) + bb*yy0(0:3) + cc*zz0(0:3)) / t1
!投影する
zz0(0:3)=zz0(0:3)+rotray(2)*t(0:3)
do j4=0, 3
if(zz0(j4) .le. cutline ) then
!カットラインより後ろにある
if( dmax1(zz0(0), zz0(1), zz0(2), zz0(3)) .le. cutline ) goto 510 !全部後ろにある場合、何もしない
!どれかは前にあり、どれかは後ろにあるので、影の３次元座標を切り取る
zz0(4:7)=zz0(0:3)
xx0(0:3)=xx0(0:3)+rotray(0)*t(0:3)
yy0(0:3)=yy0(0:3)+rotray(1)*t(0:3)
xx0(4:7)=xx0(0:3)
yy0(4:7)=yy0(0:3)
if(zz0(0) .le. cutline ) then
j9=0
do while( zz0(j9) .le. cutline )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
else
j5=0
do while( zz0(j5) .gt. cutline )
j5=j5+1
end do
j9=j5
do while( zz0(j9) .le. cutline )
j9=j9+1
end do
!切り取る線 zz0(j9-1) - zz0(j9)
j5=j9
do while( zz0(j5) .gt. cutline ) !前にある座標
j5=j5+1
end do
!切り取る線 zz0(j5-1) - zz0(j5)
end if
!zz0(j9-1)とzz0(j9)の切り取り点, 点zz0(j9)～zz0(j5-1), zz0(j5-1)とzz0(j5)の切り取り点

dz= zz0(j9)-zz0(j9-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 510 !zz0(j9-1)とzz0(j9)が近い場合、何もしない
dz= ( cutline-zz0(j9-1) ) / dz
xx0(j9-1)= xx0(j9-1) + ( xx0(j9)-xx0(j9-1) ) * dz
yy0(j9-1)= yy0(j9-1) + ( yy0(j9)-yy0(j9-1) ) * dz
zz0(j9-1)= 0.0_8+cutline

dz= zz0(j5)-zz0(j5-1)
if( dabs(dz) .lt. 0.0001_8 ) goto 510 !zz0(j5-1)とzz0(j5)が近い場合、何もしない
dz= ( cutline-zz0(j5-1) ) / dz
xx0(j5)= xx0(j5-1) + ( xx0(j5)-xx0(j5-1) ) * dz
yy0(j5)= yy0(j5-1) + ( yy0(j5)-yy0(j5-1) ) * dz
zz0(j5)= 0.0_8+cutline

zz0(j9-1:j5) = (screen+eyeoffset)/(zz0(j9-1:j5)+eyeoffset)
n= j5 - (j9-1) + 1
polydata(0:n-1).wx = xx0(j9-1:j5)*zz0(j9-1:j5)*zoom
polydata(0:n-1).wy = yy0(j9-1:j5)*zz0(j9-1:j5)*zoom
goto 430
end if
end do
!４点すべてカットラインより前にあった
zz0(0:3) = (screen+eyeoffset)/(zz0(0:3)+eyeoffset)
polydata(0:3).wx = (xx0(0:3) + rotray(0)*t(0:3))*zz0(0:3)*zoom
polydata(0:3).wy = (yy0(0:3) + rotray(1)*t(0:3))*zz0(0:3)*zoom
n=4

430 continue

if( isrgnsetted ) then
call drawpolygon2( n )
if( afterpolydots .ge. 3 ) shadowcount=shadowcount+1 !３角形以上のとき、採用された
end if

end if
510 continue
end do
500 continue
end if
end do
!影の原因になる可能性のあるすべての立方体を調べた
400 continue

end if

!画像をキャプチャーする
j2=beforeptrcapture
do while (j2 .lt. ptrcapture)
rgnxm=capture(j2+0)
rgnym=capture(j2+1)
dx=capture(j2+2)
dy=capture(j2+3)
j2=j2+4
c=dx*dy
do j3=rgnym, rgnym+dy-1 !ポリゴンをビットマスク処理したものをキャプチャーする
capture(c+j2:c+j2+dx-1)= ( vram2(rgnxm:rgnxm+dx-1,j3) .and. .not. capture(j2:j2+dx-1) )
j2=j2+dx
end do
j2=j2+c
end do

! 1グループの処理終わり
end do

!ソートされた順に実際に表示する
do j1=0, ppsort-1
j2=groupforserial(sortc(j1))
if( j2 .ne. -1 ) then
rgnxm=capture(j2+0)
rgnym=capture(j2+1)
dx=capture(j2+2)
dy=capture(j2+3)
j2=j2+4
c=dx*dy
do j3=rgnym, rgnym+dy-1
vram(rgnxm:rgnxm+dx-1,j3)= ( ( vram(rgnxm:rgnxm+dx-1,j3) .and. capture(j2:j2+dx-1) ) .or. capture(c+j2:c+j2+dx-1) )
j2=j2+dx
end do
end if
end do

end if

end if

999 continue

if(debugmode .and. shadowmode) then
capturepercent=int2(ptrcapture*100 / maxcapture)
cutanglespercent=int2(endofcutangles*100 / maxcutangles)
end if

end subroutine displaycube

integer(2) function cmp_function_polygon(a1, a2)
use common_variables

implicit none
integer(2) a1, a2
real(8) b1, b2

b1=distance(a1)
b2=distance(a2)

if( b1 .gt. b2 ) then
cmp_function_polygon=int2(-1)
return
end if

if( b1 .lt. b2 ) then
cmp_function_polygon=int2(1)
return
end if

cmp_function_polygon=int2(0)
return

end function cmp_function_polygon

!平面(sur_x0, sur_y0, sur_z0)-(sur_x1, sur_y1, sur_z1)-(sur_x2, sur_y2, sur_z2)の式
! a*x + b*y + c*z + d == 0を求める
!結果は、sur_a, sur_b, sur_c, sur_d に入る
subroutine getabcd()
use common_variables
implicit none

real(8) a(0:2, 0:5)
real(8) tmp, m, b
integer(2) j, k, l, kp

a(0,0)=sur_x0
a(0,1)=sur_x1
a(0,2)=sur_x2
a(1,0)=sur_y0
a(1,1)=sur_y1
a(1,2)=sur_y2
a(2,0)=sur_z0
a(2,1)=sur_z1
a(2,2)=sur_z2

do j=0, 2
do k=0, 2
if(j .eq. k) then
a(j,3+k)=-1.0_8
else
a(j,3+k)=0.0_8
end if
end do
end do

do j=0, 2
do k=0, 1
a(j,k)=a(j,k)-a(j,2)
end do
end do

do j=0, 1
m=dabs(a(j,j))
kp=j
do k=j+1, 2
if( dabs(a(k,j)) .gt. m ) then
m=dabs(a(k,j))
kp=k
end if
end do
if(m .eq. 0) then
result=messageboxqq('getabcdエラーです(最大ピボットが０)'C, 'Neo cubes'C, MB$OK)
goto 240
end if
if(kp .ne. j) then
do k=j, 5
tmp=a(j,k)
a(j,k)=a(kp,k)
a(kp,k)=tmp
end do
end if

do k=j+1, 2
b=a(k,j)/a(j,j)
do l=j+1, 5
a(k,l)=a(k,l)-b*a(j,l)
end do
end do

240 continue

end do

sur_a=a(2,3)
sur_b=a(2,4)
sur_c=a(2,5)
sur_d=a(2,2)

!デバッグ*******
b=dabs( sur_a*sur_x0 + sur_b*sur_y0 + sur_c*sur_z0 + sur_d )
b=b+dabs( sur_a*sur_x1 + sur_b*sur_y1 + sur_c*sur_z1 + sur_d )
b=b+dabs( sur_a*sur_x2 + sur_b*sur_y2 + sur_c*sur_z2 + sur_d )
if( b .gt. 0.0001_8 ) then
result=messageboxqq('getabcdエラーです(式が０にならない)'C, 'Neo cubes'C, MB$OK)
end if

end subroutine getabcd

subroutine drawpolygon2(d)
use common_variables
implicit none

integer(2) d
integer(2) x0, y0, x1, y1
integer(2) j4, j5
integer(2) pxx, pyy

plx(0:d-1)=polydata(0:d-1).wx
ply(0:d-1)=polydata(0:d-1).wy
spcbuf=0
dots=d
call poly2()

afterpolydots=spcbuf

if(spcbuf .lt. 3) return !３角形以下では何もしない

!点を整数化し、重複する点を除外する
!すべての点を含む矩形の範囲を知る
j4=0
x0=size_x/2 + int2(cxbuf(ptrcbuf(spcbuf-1)))
rgnxm2=x0
rgnxp2=x0
y0=size_y/2 + int2(cybuf(ptrcbuf(spcbuf-1)))
rgnym2=y0
rgnyp2=y0
do j5=0, spcbuf-1
x1=size_x/2 + int2(cxbuf(ptrcbuf(j5)))
if( x1 .lt. rgnxm2 ) rgnxm2=x1
if( x1 .gt. rgnxp2 ) rgnxp2=x1
y1=size_y/2 + int2(cybuf(ptrcbuf(j5)))
if( y1 .lt. rgnym2 ) rgnym2=y1
if( y1 .gt. rgnyp2 ) rgnyp2=y1
if( (x1 .ne. x0) .or. (y1 .ne. y0) ) then
dispdata(j4).xcoord=x1
dispdata(j4).ycoord=y1
j4=j4+1
x0=x1
y0=y1
end if
end do

afterpolydots=j4

if(j4 .lt. 3) return !３角形以下の場合は何もしない

!リージョンが明らかに重ならないときは、何もしない
if( rgnxpgroup .lt. rgnxm2 ) goto 999
if( rgnypgroup .lt. rgnym2 ) goto 999
if( rgnxp2 .lt. rgnxmgroup ) goto 999
if( rgnyp2 .lt. rgnymgroup ) goto 999

pxx=dispdata(0).xcoord
pyy=dispdata(0).ycoord

do j5=1, afterpolydots-2
trix(0)=pxx
triy(0)=pyy
trix(1)=dispdata(j5).xcoord
triy(1)=dispdata(j5).ycoord
trix(2)=dispdata(j5+1).xcoord
triy(2)=dispdata(j5+1).ycoord
call drawtriangle2()
end do

return

999 continue
afterpolydots=0

end subroutine drawpolygon2

subroutine drawtriangle2()
use common_variables
implicit none

integer(4) co
integer(2) tmp
real(4) dxybig, dxysmall
real(4) xbig, xsmall
integer(2) left, right
integer(2) j
integer(2) dybig, dysmall
real(4) d

!y座標が小さい順に並べる
if(triy(0) .gt. triy(1)) then
tmp=triy(0)
triy(0)=triy(1)
triy(1)=tmp
tmp=trix(0)
trix(0)=trix(1)
trix(1)=tmp
end if
if(triy(0) .gt. triy(2)) then
tmp=triy(0)
triy(0)=triy(2)
triy(2)=tmp
tmp=trix(0)
trix(0)=trix(2)
trix(2)=tmp
end if
if(triy(1) .gt. triy(2)) then
tmp=triy(2)
triy(2)=triy(1)
triy(1)=tmp
tmp=trix(2)
trix(2)=trix(1)
trix(1)=tmp
end if

dybig=triy(2)-triy(0)
dxybig=(trix(2)-trix(0))/(dybig*2.0_4)

!上半分
dysmall=triy(1)-triy(0)
dxysmall=(trix(1)-trix(0))/(dysmall*2.0_4)
select case (dysmall)

case(0) !横一線
left=imin0( trix(0), trix(1) )
right=imax0( trix(0), trix(1) )
vram2(left:right, triy(0))= ( vram2(left:right, triy(0)) .and. shadowcolor ) !塗る

case default !一般
xbig=trix(0)+dxybig
xsmall=trix(0)+dxysmall
!最初の横一線
left=int2( 0.5_4+ amin1( trix(0)+0.0_4, xbig, xsmall ) )
right=int2( 0.5_4+ amax1( trix(0)+0.0_4, xbig, xsmall ) )
vram2(left:right, triy(0))= ( vram2(left:right, triy(0)) .and. shadowcolor ) !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(0)+1, triy(1)-1
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
vram2(left:right, j)= ( vram2(left:right, j) .and. shadowcolor ) !塗る
xbig=xbig+dxybig*2
xsmall=xsmall+dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

end select

!下半分
dysmall=triy(2)-triy(1)
dxysmall=(trix(2)-trix(1))/(dysmall*2.0_4)
select case (dysmall)

case(0) !横一線
left=imin0( trix(1), trix(2) )
right=imax0( trix(1), trix(2) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

case default !一般
xbig=trix(2)-dxybig
xsmall=trix(2)-dxysmall
!最初の横一線
left=int2( 0.5_4+ amin1( trix(2)+0.0_4, xbig, xsmall ) )
right=int2( 0.5_4+ amax1( trix(2)+0.0_4, xbig, xsmall ) )
vram2(left:right, triy(2))= ( vram2(left:right, triy(2)) .and. shadowcolor ) !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(2)-1, triy(1)+1, -1
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
vram2(left:right, j)= ( vram2(left:right, j) .and. shadowcolor ) !塗る
xbig=xbig-dxybig*2
xsmall=xsmall-dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
vram2(left:right, triy(1))= ( vram2(left:right, triy(1)) .and. shadowcolor ) !塗る

end select

end subroutine drawtriangle2

!(polydata(0).wx,polydata(0).wy) - (polydata(1).wx,polydata(1).wy) - ... -
! (polydata(d-1).wx,polydata(d-1).wy) を色 co で vram にポリゴンを表示する
!afterpolydots=実際に描画された多角形の頂点数
subroutine drawpolygon(co, d)
use common_variables
implicit none

integer(4) co
integer(2) d
integer(2) x0, y0, x1, y1
integer(2) j4, j5
integer(2) pxx, pyy

plx(0:d-1)=polydata(0:d-1).wx
ply(0:d-1)=polydata(0:d-1).wy
spcbuf=0
dots=d
call poly()

afterpolydots=spcbuf

if(spcbuf .lt. 3) goto 999 !３角形以下では何もしない

!点を整数化し、重複する点を除外する
!すべての点を含む矩形の範囲を知る
j4=0
x0=size_x/2 + int2(cxbuf(ptrcbuf(spcbuf-1)))
rgnxm=x0
rgnxp=x0
y0=size_y/2 + int2(cybuf(ptrcbuf(spcbuf-1)))
rgnym=y0
rgnyp=y0
do j5=0, spcbuf-1
x1=size_x/2 + int2(cxbuf(ptrcbuf(j5)))
if( x1 .lt. rgnxm ) rgnxm=x1
if( x1 .gt. rgnxp ) rgnxp=x1
y1=size_y/2 + int2(cybuf(ptrcbuf(j5)))
if( y1 .lt. rgnym ) rgnym=y1
if( y1 .gt. rgnyp ) rgnyp=y1
if( (x1 .ne. x0) .or. (y1 .ne. y0) ) then
dispdata(j4).xcoord=x1
dispdata(j4).ycoord=y1
j4=j4+1
x0=x1
y0=y1
end if
end do

afterpolydots=j4

if(j4 .lt. 3) goto 999 !３角形以下の場合は何もしない

!グループ全体のリージョンを設定する
if( .not. isrgnsetted ) then
isrgnsetted=.true. !リージョンが１回以上セットされた
rgnxmgroup=rgnxm
rgnxpgroup=rgnxp
rgnymgroup=rgnym
rgnypgroup=rgnyp
else
if( rgnxm .lt. rgnxmgroup ) rgnxmgroup=rgnxm
if( rgnxp .gt. rgnxpgroup ) rgnxpgroup=rgnxp
if( rgnym .lt. rgnymgroup ) rgnymgroup=rgnym
if( rgnyp .gt. rgnypgroup ) rgnypgroup=rgnyp
end if

vram3(rgnxm:rgnxp, rgnym:rgnyp)=16#ffffff

pxx=dispdata(0).xcoord
pyy=dispdata(0).ycoord

do j5=1, afterpolydots-2
trix(0)=pxx
triy(0)=pyy
trix(1)=dispdata(j5).xcoord
triy(1)=dispdata(j5).ycoord
trix(2)=dispdata(j5+1).xcoord
triy(2)=dispdata(j5+1).ycoord
call drawtriangle(co)
end do

999 continue

end subroutine drawpolygon

!整数座標( trix(0), triy(0) ) - ( trix(1), triy(1) ) - ( trix(2), triy(2) )
!の３角形を色 co で vram に表示する
subroutine drawtriangle(co)
use common_variables
implicit none

integer(4) co
integer(2) tmp
real(4) dxybig, dxysmall
real(4) xbig, xsmall
integer(2) left, right
integer(2) j
integer(2) dybig, dysmall
real(4) d

!y座標が小さい順に並べる
if(triy(0) .gt. triy(1)) then
tmp=triy(0)
triy(0)=triy(1)
triy(1)=tmp
tmp=trix(0)
trix(0)=trix(1)
trix(1)=tmp
end if
if(triy(0) .gt. triy(2)) then
tmp=triy(0)
triy(0)=triy(2)
triy(2)=tmp
tmp=trix(0)
trix(0)=trix(2)
trix(2)=tmp
end if
if(triy(1) .gt. triy(2)) then
tmp=triy(2)
triy(2)=triy(1)
triy(1)=tmp
tmp=trix(2)
trix(2)=trix(1)
trix(1)=tmp
end if

dybig=triy(2)-triy(0)
dxybig=(trix(2)-trix(0))/(dybig*2.0_4)

!上半分
dysmall=triy(1)-triy(0)
dxysmall=(trix(1)-trix(0))/(dysmall*2.0_4)
select case (dysmall)

case(0) !横一線
left=imin0( trix(0), trix(1) )
right=imax0( trix(0), trix(1) )
vram2(left:right, triy(0))= co !塗る
vram3(left:right, triy(0))= 0 !塗る

case default !一般
xbig=trix(0)+dxybig
xsmall=trix(0)+dxysmall
!最初の横一線
left=int2( 0.5_4+ amin1( trix(0)+0.0_4, xbig, xsmall ) )
right=int2( 0.5_4+ amax1( trix(0)+0.0_4, xbig, xsmall ) )
vram2(left:right, triy(0))= co !塗る
vram3(left:right, triy(0))= 0 !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(0)+1, triy(1)-1
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig+dxybig*2, xsmall+dxysmall*2 ) )
vram2(left:right, j)= co !塗る
vram3(left:right, j)= 0 !塗る
xbig=xbig+dxybig*2
xsmall=xsmall+dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig+d, xsmall+dxysmall ) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

end select

!下半分
dysmall=triy(2)-triy(1)
dxysmall=(trix(2)-trix(1))/(dysmall*2.0_4)
select case (dysmall)

case(0) !横一線
left=imin0( trix(1), trix(2) )
right=imax0( trix(1), trix(2) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

case default !一般
xbig=trix(2)-dxybig
xsmall=trix(2)-dxysmall
!最初の横一線
left=int2( 0.5_4+ amin1( trix(2)+0.0_4, xbig, xsmall ) )
right=int2( 0.5_4+ amax1( trix(2)+0.0_4, xbig, xsmall ) )
vram2(left:right, triy(2))= co !塗る
vram3(left:right, triy(2))= 0 !塗る
!中の線
if(dysmall .ge. 2) then
do j=triy(2)-1, triy(1)+1, -1
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig-dxybig*2, xsmall-dxysmall*2 ) )
vram2(left:right, j)= co !塗る
vram3(left:right, j)= 0 !塗る
xbig=xbig-dxybig*2
xsmall=xsmall-dxysmall*2
end do
end if
!最後の横一線
if(dybig .eq. dysmall) then
d=dxybig
else
d=dxybig*2
end if
left=int2( 0.5_4+ amin1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
right=int2( 0.5_4+ amax1( xbig, xsmall, xbig-d, xsmall-dxysmall ) )
vram2(left:right, triy(1))= co !塗る
vram3(left:right, triy(1))= 0 !塗る

end select

end subroutine drawtriangle

!(plx(0), ply(0)) - (ply(1), ply(1)) - ... - (plx(dots-1), ply(dots-1))
! のポリゴンと四角リージョンの論理積をとる
!cxbuf(spcbuf++), cybuf(spcbuf++) に atan2 でソートされた頂点が入る
!ソートされた結果はポインタの配列 ptrcbuf(0:spcbuf-1) に入っている
subroutine poly()
use common_variables
implicit none
integer(2), external :: cmp_function_poly
real(8) poly_a(0:dots-1), poly_b(0:dots-1), poly_c(0:dots-1)
real(8) poly_sign(0:dots-1)
real(8) cx, cy, sign
integer(2) j, k

!四角リージョンの内側にあるポリゴンの頂点は採用する
k=0
do j=0, dots-1
if( plx(j) .ge. 0.0_8-size_x/2 ) then
if( plx(j) .le. size_x/2-1.0_8 ) then
if( ply(j) .ge. 0.0_8-size_y/2 ) then
if( ply(j) .le. size_y/2-1.0_8 ) then
cxbuf(spcbuf)=plx(j)
cybuf(spcbuf)=ply(j)
spcbuf=spcbuf+1
k=k+1
end if
end if
end if
end if
end do

if(k .eq. dots) goto 100 !すべての点が四角リージョンの内側にある場合、すぐにソート処理へ行く

cx=0.0_8
cy=0.0_8
do j=0, dots-1
cx=cx+plx(j)
cy=cy+ply(j)
end do
cx=cx/dots
cy=cy/dots

do j=0, dots-2
line_x0=plx(j)
line_y0=ply(j)
line_x1=plx(j+1)
line_y1=ply(j+1)
call getabc()
poly_a(j)=line_a
poly_b(j)=line_b
poly_c(j)=line_c
end do
line_x0=plx(dots-1)
line_y0=ply(dots-1)
line_x1=plx(0)
line_y1=ply(0)
call getabc()
poly_a(dots-1)=line_a
poly_b(dots-1)=line_b
poly_c(dots-1)=line_c

poly_sign(0:dots-1)= poly_a(0:dots-1)*cx + poly_b(0:dots-1)*cy + poly_c(0:dots-1)

!ポリゴンの内側にある四角リージョンの頂点は採用する
do j=0, 3
do k=0, dots-1
sign= poly_a(k)*rctx(j) + poly_b(k)*rcty(j) + poly_c(k)
if( poly_sign(k)*sign .lt. 0.0_8 ) goto 200 !異符号は採用しない
end do
cxbuf(spcbuf)=rctx(j)
cybuf(spcbuf)=rcty(j)
spcbuf=spcbuf+1
200 continue
end do

!ポリゴンの線と四角リージョンの線の交点が有効なものは採用する
do j=0, dots-2
crx(0)=plx(j)
cry(0)=ply(j)
crx(1)=plx(j+1)
cry(1)=ply(j+1)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()
end do
crx(0)=plx(dots-1)
cry(0)=ply(dots-1)
crx(1)=plx(0)
cry(1)=ply(0)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()

100 continue

!ソートする
if(spcbuf .gt. 0) then
do j=0, spcbuf-1
ptrcbuf(j)=j
end do
end if

if(spcbuf .gt. 3) then

cx=0.0_8
cy=0.0_8
do j=0, spcbuf-1
cx=cx+cxbuf(j)
cy=cy+cybuf(j)
end do
cx=cx/spcbuf
cy=cy/spcbuf
do j=0, spcbuf-1
degree_cbuf(j)= datan2( cybuf(j)-cy, cxbuf(j)-cx )
end do
CALL qsort(ptrcbuf, spcbuf, 2, cmp_function_poly)

end if

end subroutine poly

integer(2) function cmp_function_poly(a1, a2)
use common_variables

implicit none
integer(2) a1, a2
real(8) b1, b2

b1=degree_cbuf(a1)
b2=degree_cbuf(a2)

if( b1 .gt. b2 ) then
cmp_function_poly=int2(-1)
return
end if

if( b1 .lt. b2 ) then
cmp_function_poly=int2(1)
return
end if

cmp_function_poly=int2(0)
return

end function cmp_function_poly

subroutine poly2()
use common_variables
implicit none
integer(2), external :: cmp_function_poly
real(8) poly_a(0:dots-1), poly_b(0:dots-1), poly_c(0:dots-1)
real(8) poly_sign(0:dots-1)
real(8) cx, cy, sign
integer(2) j, k
real(8) rgnxmr, rgnxpr, rgnymr, rgnypr, x, y

!四角リージョンの内側にあるポリゴンの頂点は採用する
rgnxmr=plx(0)
rgnxpr=rgnxmr
rgnymr=ply(0)
rgnypr=rgnymr
k=0
do j=0, dots-1
x=plx(j)
if( x .lt. rgnxmr ) rgnxmr=x
if( x .gt. rgnxpr ) rgnxpr=x
y=ply(j)
if( y .lt. rgnymr ) rgnymr=y
if( y .gt. rgnypr ) rgnypr=y

if( x .ge. 0.0_8-size_x/2 ) then
if( x .le. size_x/2-1.0_8 ) then
if( y .ge. 0.0_8-size_y/2 ) then
if( y .le. size_y/2-1.0_8 ) then
cxbuf(spcbuf)=x
cybuf(spcbuf)=y
spcbuf=spcbuf+1
k=k+1
end if
end if
end if
end if
end do

!リージョンが明らかに重ならないときは、何もしない
if( 0.0_8+rgnxpgroup .lt. size_x/2+rgnxmr ) goto 999
if( 0.0_8+rgnypgroup .lt. size_y/2+rgnymr ) goto 999
if( size_x/2+rgnxpr .lt. 0.0_8+rgnxmgroup ) goto 999
if( size_y/2+rgnypr .lt. 0.0_8+rgnymgroup ) goto 999

if(k .eq. dots) goto 100 !すべての点が四角リージョンの内側にある場合、すぐにソート処理へ行く

cx=0.0_8
cy=0.0_8
do j=0, dots-1
cx=cx+plx(j)
cy=cy+ply(j)
end do
cx=cx/dots
cy=cy/dots

do j=0, dots-2
line_x0=plx(j)
line_y0=ply(j)
line_x1=plx(j+1)
line_y1=ply(j+1)
call getabc()
poly_a(j)=line_a
poly_b(j)=line_b
poly_c(j)=line_c
end do
line_x0=plx(dots-1)
line_y0=ply(dots-1)
line_x1=plx(0)
line_y1=ply(0)
call getabc()
poly_a(dots-1)=line_a
poly_b(dots-1)=line_b
poly_c(dots-1)=line_c

poly_sign(0:dots-1)= poly_a(0:dots-1)*cx + poly_b(0:dots-1)*cy + poly_c(0:dots-1)

!ポリゴンの内側にある四角リージョンの頂点は採用する
do j=0, 3
do k=0, dots-1
sign= poly_a(k)*rctx(j) + poly_b(k)*rcty(j) + poly_c(k)
if( poly_sign(k)*sign .lt. 0.0_8 ) goto 200 !異符号は採用しない
end do
cxbuf(spcbuf)=rctx(j)
cybuf(spcbuf)=rcty(j)
spcbuf=spcbuf+1
200 continue
end do

!ポリゴンの線と四角リージョンの線の交点が有効なものは採用する
do j=0, dots-2
crx(0)=plx(j)
cry(0)=ply(j)
crx(1)=plx(j+1)
cry(1)=ply(j+1)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()
end do
crx(0)=plx(dots-1)
cry(0)=ply(dots-1)
crx(1)=plx(0)
cry(1)=ply(0)
do k=0, 2
crx(2)=rctx(k)
cry(2)=rcty(k)
crx(3)=rctx(k+1)
cry(3)=rcty(k+1)
call cross()
end do
crx(2)=rctx(3)
cry(2)=rcty(3)
crx(3)=rctx(0)
cry(3)=rcty(0)
call cross()

100 continue

!ソートする
if(spcbuf .gt. 0) then
do j=0, spcbuf-1
ptrcbuf(j)=j
end do
end if

if(spcbuf .gt. 3) then

cx=0.0_8
cy=0.0_8
do j=0, spcbuf-1
cx=cx+cxbuf(j)
cy=cy+cybuf(j)
end do
cx=cx/spcbuf
cy=cy/spcbuf
do j=0, spcbuf-1
degree_cbuf(j)= datan2( cybuf(j)-cy, cxbuf(j)-cx )
end do
CALL qsort(ptrcbuf, spcbuf, 2, cmp_function_poly)

end if

return

999 continue
spcbuf=0

end subroutine poly2

!(0, 0) - (size_x-1, size_y-1) の四角リージョンを設定する（自作シェーダー表示用）
!最初に１度だけ呼ばれる
!rctrgn_a(0~3), rctrgn_b(0~3), rctrgn_c(0~3) に各直線のa,b,cが入る
!rctrgn_sign(0~3) にポリゴンの真ん中の点に対する符号が入る
subroutine setrctrgn()
use common_variables
implicit none
real(8) cx, cy
integer(2) j

rctx(0)=0.0_8-size_x/2
rcty(0)=0.0_8-size_y/2
rctx(1)=0.0_8+size_x/2-1
rcty(1)=0.0_8-size_y/2
rctx(2)=0.0_8+size_x/2-1
rcty(2)=0.0_8+size_y/2-1
rctx(3)=0.0_8-size_x/2
rcty(3)=0.0_8+size_y/2-1

cx=0.0_8
cy=0.0_8
do j=0, 3
cx=cx+rctx(j)
cy=cy+rcty(j)
end do
cx=cx/4
cy=cy/4

do j=0, 2
line_x0=rctx(j)
line_y0=rcty(j)
line_x1=rctx(j+1)
line_y1=rcty(j+1)
call getabc()
rctrgn_a(j)=line_a
rctrgn_b(j)=line_b
rctrgn_c(j)=line_c
end do
line_x0=rctx(3)
line_y0=rcty(3)
line_x1=rctx(0)
line_y1=rcty(0)
call getabc()
rctrgn_a(3)=line_a
rctrgn_b(3)=line_b
rctrgn_c(3)=line_c

rctrgn_sign(0:3)= rctrgn_a(0:3)*cx + rctrgn_b(0:3)*cy + rctrgn_c(0:3)

end subroutine setrctrgn

!直線(line_x0, line_y0)-(line_x1, line_y1)の式 a*x + b*y + c == 0を求める
!結果は、line_a, line_b, line_c に入る
subroutine getabc()
use common_variables
implicit none

!dabs(x1-x0) と dabs(y1-y0) のスパンの大きいほうを選ぶ
!
!y1-y0が大きいとき、次の方程式を解く
!(x0 + b*y0 + c == 0) && (x1 + b*y1 + c == 0)
!a=1, b=(x0 - x1)/(y1-y0), c=(x1 y0 - x0 y1)/(y1-y0)
!
!x1-x0が大きいとき、次の方程式を解く
!(a*x0 + y0 + c == 0) && (a*x1 + y1 + c == 0)
!a=(y0 - y1)/(x1-x0), b=1, c=(x0 y1 -x1 y0)/(x1-x0)

real(8) x1mx0, y1my0
real(8) x1y0mx0y1

x1mx0=line_x1-line_x0
y1my0=line_y1-line_y0
x1y0mx0y1=line_x1*line_y0-line_x0*line_y1

if( dabs(x1mx0) .gt. dabs(y1my0) ) then
!x1-x0が大きいとき
if(x1mx0 .eq. 0.0_8) then
!点になっている
line_a=1.0_8
line_b=0.0_8
line_c=-line_x0
else
!a=(y0 - y1)/(x1-x0), b=1, c=(x0 y1 -x1 y0)/(x1-x0)
line_a=-y1my0/x1mx0
line_b=1.0_8
line_c=-x1y0mx0y1/x1mx0
end if
else
!y1-y0が大きいとき
if(y1my0 .eq. 0.0_8) then
!点になっている
line_a=0.0_8
line_b=1.0_8
line_c=-line_y0
else
!a=1, b=(x0 - x1)/(y1-y0), c=(x1 y0 - x0 y1)/(y1-y0)
line_a=1.0_8
line_b=-x1mx0/y1my0
line_c=x1y0mx0y1/y1my0
endif
end if

end subroutine getabc

!直線(crx(0), cry(0))-(crx(1), cry(1)) と (crx(2), cry(2))-(crx(3), cry(3))の交点を計算する
!結果は、cxbuf(spcbuf++), cybuf(spcbuf++)に入る
subroutine cross()
use common_variables
implicit none

!次の連立方程式を s 及び t について解く
!
!x=x0+( x1-x0 )*s
!y=y0+( y1-y0 )*s
!
!x=x2+( x3-x2 )*t
!y=y2+( y3-y2 )*t
!
!s=x3 (y0 - y2) + x0 (y2 - y3) + x2 (-y0 + y3) /
! -(x2 - x3) (y0 - y1) + (x0 - x1) (y2 - y3)
!
!t=x2 (y0 - y1) + x0 (y1 - y2) + x1 (-y0 + y2) /
! (x2 - x3) (y0 - y1) - (x0 - x1) (y2 - y3)
!
!0<=s<=1 .and. 0<=t<=1 ならば交点が存在する

real(8) x0,y0, x1,y1, x2,y2, x3,y3
real(8) s, t, u, us, ut
real(8) y0y1, y0y2, y2y3
real(8) tmp, abs01, abs23, xm, ym, xp, yp
logical(4) reverse

x0=crx(0)
y0=cry(0)
x1=crx(1)
y1=cry(1)
x2=crx(2)
y2=cry(2)
x3=crx(3)
y3=cry(3)

y0y1=y0-y1
y0y2=y0-y2
y2y3=y2-y3

u= (x2 - x3) * y0y1 - (x0 - x1) * y2y3
us= x3 * y0y2 + x0 * y2y3 + x2 * (y3 - y0)
ut= x2 * y0y1 + x0 * (y1 - y2) - x1 * y0y2

if( u .ne. 0.0_8 ) then

s= us / -u
t= ut / u

if( s .ge. 0.0_8 ) then
if( s .le. 1.0_8 ) then
if( t .ge. 0.0_8 ) then
if( t .le. 1.0_8 ) then
cxbuf(spcbuf)=x0+(x1-x0)*s
cybuf(spcbuf)=y0+(y1-y0)*s
spcbuf=spcbuf+1
end if
end if
end if
end if

else

!u==0のとき
if( (us .eq. 0.0_8) .and. (ut .eq. 0.0_8) ) then
!２つの線が１直線上に重なっている場合を救済する
abs01=(x1-x0)**2+(y1-y0)**2
abs23=(x3-x2)**2+(y3-y2)**2
!スパンの長い方をx0,y0,x1,y1とする
if( abs01 .lt. abs23 ) then
abs01=abs23
tmp=x0
x0=x2
x2=tmp
tmp=y0
y0=y2
y2=tmp
tmp=x1
x1=x3
x3=tmp
tmp=y1
y1=y3
y3=tmp
end if
if( abs01 .eq. 0.0_8 ) then
!点になっている
if( (x0 .eq. x2) .and. (y0 .eq. y2) ) then
cxbuf(spcbuf)=x0
cybuf(spcbuf)=y0
spcbuf=spcbuf+1
end if
else
!スパンの長い方をxとする
if( dabs(x1-x0) .lt. dabs(y1-y0) ) then
tmp=x0
x0=y0
y0=tmp
tmp=x1
x1=y1
y1=tmp
tmp=x2
x2=y2
y2=tmp
tmp=x3
x3=y3
y3=tmp
reverse=.true.
else
reverse=.false.
end if
!直線で比較する
if( x0 .gt. x1 ) then
tmp=x0
x0=x1
x1=tmp
tmp=y0
y0=y1
y1=tmp
end if
if( x2 .gt. x3 ) then
tmp=x2
x2=x3
x3=tmp
tmp=y2
y2=y3
y3=tmp
end if
!マイナス方向をカット
if( x0 .lt. x2 ) then
xm=x2
ym=y2
else
xm=x0
ym=y0
end if
!プラス方向をカット
if( x1 .gt. x3 ) then
xp=x3
yp=y3
else
xp=x1
yp=y1
end if
!---> xm,ym ---> xp,yp ---> の場合、救済可能
if( xm .le. xp ) then
!(xm,ym)および(xp,yp)を登録
if(reverse) then
tmp=xm
xm=ym
ym=tmp
tmp=xp
xp=yp
yp=tmp
end if
cxbuf(spcbuf)=xm
cybuf(spcbuf)=ym
spcbuf=spcbuf+1
cxbuf(spcbuf)=xp
cybuf(spcbuf)=yp
spcbuf=spcbuf+1
end if
end if
end if

end if

end subroutine cross

!ファイルをロードする
subroutine loadfile()
use common_variables
USE IFPORT
implicit none

CHARACTER(256) pathname
INTEGER(4) pathlen
integer(4) io
integer(4) j

result = SETACTIVEQQ (childwindow1)

CALL SETTEXTPOSITION (INT2(1), INT2(1), curpos)
write(childwindow1,*) '＊＊＊　ロード　＊＊＊'
write(childwindow1,*) 'ファイルをロードします'
write(childwindow1,*) '1～9, a～z のばんごうを押してください（0=やめる）'

call revewfile()

100 continue
call selectfilenumber()
if(filenumber .eq. -1) goto 100

if(filenumber .ne. 0) then

pathlen = FINDFILEQQ(filenamelist(filenumber), "", pathname)
if (pathlen .eq. 0) then
result=messageboxqq('ファイルがありません'C, 'Neo cubes'C, MB$OK)
goto 999
end if

open(fileequipnumber, file=filenamelist(filenumber), iostat=io)
if (io .ne. 0) then
result=messageboxqq('オープンできません'C, 'Neo cubes'C, MB$OK)
goto 999
end if

!読み込む
read(fileequipnumber, *, end=888) ptrcube !立方体の位置データ数

read(fileequipnumber, *, end=888) eye(0), eye(1), eye(2) !現在位置

read(fileequipnumber, *, end=888) matrix0(0,0) !視点方向
read(fileequipnumber, *, end=888) matrix0(1,0) !視点方向
read(fileequipnumber, *, end=888) matrix0(2,0) !視点方向
read(fileequipnumber, *, end=888) matrix0(0,1) !視点方向
read(fileequipnumber, *, end=888) matrix0(1,1) !視点方向
read(fileequipnumber, *, end=888) matrix0(2,1) !視点方向
read(fileequipnumber, *, end=888) matrix0(0,2) !視点方向
read(fileequipnumber, *, end=888) matrix0(1,2) !視点方向
read(fileequipnumber, *, end=888) matrix0(2,2) !視点方向

read(fileequipnumber, *, end=888) ray(0) !光源方向
read(fileequipnumber, *, end=888) ray(1) !光源方向
read(fileequipnumber, *, end=888) ray(2) !光源方向

read(fileequipnumber, *, end=888) shadowmode !影モード

if(ptrcube .gt. 0) then
do j=0, ptrcube-1
!cube(xyz, 通し番号)
read(fileequipnumber, *, end=888) cube(0,j), cube(1,j), cube(2,j)
end do
end if

888 continue

close(fileequipnumber, iostat=io)
if (io .ne. 0) then
result=messageboxqq('クローズできません'C, 'Neo cubes'C, MB$OK)
goto 999
end if

result=messageboxqq('読み込みました'C, 'Neo cubes'C, MB$OK)

999 continue

end if

end subroutine loadfile

!ファイルをセーブする
subroutine savefile()
use common_variables
USE IFPORT
implicit none

integer(4) io
CHARACTER(256) pathname
INTEGER(4) pathlen
integer(4) j

result = SETACTIVEQQ (childwindow1)

CALL SETTEXTPOSITION (INT2(1), INT2(1), curpos)
write(childwindow1,*) '＊＊＊　セーブ　＊＊＊'
write(childwindow1,*) 'ファイルをセーブします'
write(childwindow1,*) '1～9, a～z のばんごうを押してください（0=やめる）'

call revewfile()

100 continue
call selectfilenumber()
if(filenumber .eq. -1) goto 100

if(filenumber .ne. 0) then

pathlen = FINDFILEQQ(filenamelist(filenumber), "", pathname)
if (pathlen .ne. 0) then
result=messageboxqq('ファイル上書きしますか？'C, 'Neo cubes'C, MB$OKCANCEL)
if (result .eq. MB$IDCANCEL) goto 999
end if

open(fileequipnumber, file=filenamelist(filenumber), iostat=io)
if (io .ne. 0) then
result=messageboxqq('オープンできません'C, 'Neo cubes'C, MB$OK)
goto 999
end if

!書き込む
write(fileequipnumber, *, iostat=io) ptrcube !立方体の位置データ数

write(fileequipnumber, *, iostat=io) eye(0), eye(1), eye(2) !現在位置

write(fileequipnumber, *, iostat=io) matrix0(0,0) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(1,0) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(2,0) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(0,1) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(1,1) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(2,1) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(0,2) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(1,2) !視点方向
write(fileequipnumber, *, iostat=io) matrix0(2,2) !視点方向

write(fileequipnumber, *, iostat=io) ray(0) !光源方向
write(fileequipnumber, *, iostat=io) ray(1) !光源方向
write(fileequipnumber, *, iostat=io) ray(2) !光源方向

write(fileequipnumber, *, iostat=io) shadowmode !影モード

if(ptrcube .gt. 0) then
do j=0, ptrcube-1
!cube(xyz, 通し番号)
write(fileequipnumber, *, iostat=io) cube(0,j), cube(1,j), cube(2,j)
end do
end if

close(fileequipnumber, iostat=io)
if (io .ne. 0) then
result=messageboxqq('クローズできません'C, 'Neo cubes'C, MB$OK)
goto 999
end if

result=messageboxqq('書き込みました'C, 'Neo cubes'C, MB$OK)

999 continue

end if

end subroutine savefile

!ファイルの有無を表示する
subroutine revewfile()
use common_variables
USE IFPORT
implicit none

CHARACTER(256) pathname
INTEGER(4) pathlen
integer(2) j
character(1) i

do j=1, 35
if( (j .ge. 1) .and. (j .le. 9) ) then
i=char(48+j)
end if
if( (j .ge. 10) .and. (j .le. 35) ) then
i=char(97+j-10)
end if
pathlen = FINDFILEQQ(filenamelist(j), "", pathname)
if (pathlen .ne. 0) then
write(childwindow1,*) 'ファイル番号', j, '(', i, ')：　満'
else
write(childwindow1,*) 'ファイル番号', j, '(', i, ')：　空'
end if
end do

end subroutine revewfile

!ファイル番号を入力する
subroutine selectfilenumber()
use common_variables
implicit none

filenumber=-1
key=incharqq()
select case (key)
case(ichar('1'))
filenumber=1
case(ichar('2'))
filenumber=2
case(ichar('3'))
filenumber=3
case(ichar('4'))
filenumber=4
case(ichar('5'))
filenumber=5
case(ichar('6'))
filenumber=6
case(ichar('7'))
filenumber=7
case(ichar('8'))
filenumber=8
case(ichar('9'))
filenumber=9
case(ichar('0'))
filenumber=0
end select

if( (key .ge. ichar('a')) .and. (key .le. ichar('z')) ) then
filenumber=key-ichar('a')+10
end if

end subroutine selectfilenumber