（その４）では，任意の超平面

　　α1Ｘ1＋α2Ｘ2＝０

と超楕円体

　　ｘ’Ｈｘ＝１

の交わりをｍ次元空間へ投影すると，ｍ次元楕円体

　　Ｘ1’［Ａ−Ｂα2~α1−（α2~α1）’Ｃ＋（α2~α1）’Ｄα2~α1］Ｘ1＝１

が得られることがわかりましたが，ｍ次元空間に投影するとはいっても，最終的に表示するのはモニタ画面上ですから平面です．したがって，（ｉ，ｊ）平面上にｍ次元楕円体を描く手順が必要とされます．

　

　その際，

　　(1)行列：Ｍ＝［Ａ−Ｂα2~α1−（α2~α1）’Ｃ＋（α2~α1）’Ｄα2~α1］

を求める

　　(2)Ｍの小行列Ｍijの固有値・固有ベクトルを求める

　　(3)軸の長さを１／√λとする

とすると，ｍ次元楕円体の（ｉ，ｊ）平面での切り口が描かれてしまいます．

　

　問題点は何かを意識していないと肝心なところがわからないので，くどくど説明しますが，いま求められているのは，ｍ次元楕円の２次元平面での切り口ではなく，ｍ次元楕円体を２次元平面に正射影することです．

　

　その計算手順は，以下のようになります．

　　(1)行列：Ｍ＝［Ａ−Ｂα2~α1−（α2~α1）’Ｃ＋（α2~α1）’Ｄα2~α1］を求める

　　(2)Ｍの逆行列Ｗを求める

　　(3)逆行列の小行列Ｗijの固有値・固有ベクトルを求める

　　(4)軸の長さを√λとする

　

　なお，３次元空間に投影する際は，

　　(3)逆行列の小行列Ｗijkの固有値・固有ベクトルを求める

ことになります．

　

　ここで，(1)の前に

　　(a)超平面のパラメータｍとα＝［α1，α2］を定める

　　(b)ムーア・ペンローズ逆行列

　　α2~＝α2’α2（α2’α2α2’α2）^(-1)α2’

あるいは

　　α2~＝α2’（α2α2’）^(-1)α2（α2’α2）^(-1)α2’

を求める

　　(c)α2~α1を求める

という手順のサブルーチンを組めば，実用的なプログラムを作成することが可能になるはずです．早速作成してみましょう．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【２】プログラム

　

1000 '

1010 ' **** robot-arm control ****

1020 '　　　　　　　　　　　　2002/01/18 (C) ｻﾄｳ ｲｸﾛｳ

1030 SCREEN 3,0:CONSOLE ,,0,1

1040 CLS 3:WIDTH 80,25:COLOR 6

1050 GOSUB *INITIALIZE

1060 '

1070　GOSUB *HYPERPLANE

1080 'GOSUB *HYPER.CONTOUR

1090 'GOSUB *HYPER.SECTION

1100　GOSUB *ELLIPSOID:GOTO 1230

1110 '

1120 STAGE$="1":' [contour]

1130 'GOSUB *HYPERPLANE

1140　GOSUB *HYPER.CONTOUR

1150 'GOSUB *ELLIPSOID.WF

1160　GOSUB *ELLIPSOID3

1170 '

1180 STAGE$="2":' [section]

1190 'GOSUB *HYPERPLANE

1200　GOSUB *HYPER.SECTION

1210 'GOSUB *ELLIPSOID.WF

1220　GOSUB *ELLIPSOID3

1230 '

1240 CLS 3

1250 END

1260 '

1270 ' *** 初期設定 ***

1280 '

1290 *INITIALIZE:

1300 RX=.7:RY=.8:BOX=1:KIND=1:POWER=0

1310 IX=.7:IY=.8:MARK=1:CLR=2:DRAW.COLOR=4

1320 TXT.COLOR=6:GRP.COLOR=5:REF.COLOR=6

1330 SCALE.KIND=1:PLOT=1

1340 PI=3.14159

1350 JX=.709:JY=1

1360 RX=IX*JX:RY=IY*JY

1370 '

1380 M=6:M0=M

1390 N1=7

1400 '

1410 MMAX=10:NMAX=10

1420 DIM V(MMAX)

1430 DIM JACOBI(MMAX,NMAX)

1440 DIM AZ(MMAX,MMAX),BZ(MMAX,MMAX)

1450 DIM HESSE(MMAX,MMAX)

1460 DIM COV(MMAX,MMAX)

1470 DIM VO(MMAX),VL(MMAX)

1480 DIM VD(MMAX),VS(MMAX)

1490 DIM UNIT$(MMAX)

1500 '

1510 DIM EDGE(MMAX)

1520 DIM S(MMAX),S0(MMAX)

1530 DIM S1(MMAX),S2(MMAX)

1540 DIM SS(MMAX*2,MMAX)

1550 DIM GXS(MMAX*2),GYS(MMAX*2)

1560 DIM GXE(MMAX*2),GYE(MMAX*2)

1570 DIM GZ(MMAX*2)

1580 '

1590 DIM XX1(MMAX,MMAX),YY1(MMAX,MMAX)

1600 DIM H1(MMAX,MMAX),H2(MMAX,MMAX)

1610 DIM H3(MMAX,MMAX)

1620 DIM C1(MMAX,MMAX),C2(MMAX,MMAX)

1630 DIM C3(MMAX,MMAX),C4(MMAX,MMAX)

1640 DIM C5(MMAX,MMAX),C6(MMAX,MMAX)

1650 DIM C7(MMAX,MMAX),C8(MMAX,MMAX)

1660 DIM ALPHA(MMAX),ALPHA1(MMAX),ALPHA2(MMAX)

1670 DIM ALPHA3(MMAX,MMAX),ALPHA4(MMAX,MMAX)

1680 '

1690 FOR I=1 TO 3

1700　 VO(I)=0:VL(I)=.5

1710　'VO(I)=0:VL(I)=1

1720　 VD(I)=5:VS(I)=1

1730　 UNIT$(I)="m/s"

1740 NEXT I

1750 FOR I=4 TO 6

1760　 VO(I)=0:VL(I)=2

1770　'VO(I)=0:VL(I)=4

1780　 VD(I)=4:VS(I)=1

1790　 UNIT$(I)="rad/s"

1800 NEXT I

1810 RETURN

1820 '

1830 ' *** 任意の超平面による切断面 ***

1840 '

1850 *HYPERPLANE:' [arbitary section]

1860 ALPHA(1)=1

1870 ALPHA(2)=1

1880 ALPHA(3)=1

1890 ALPHA(4)=1

1900 ALPHA(5)=1

1910 ALPHA(6)=1

1920 '

1930 S=0

1940 FOR I=1 TO M

1950　 S=S+ALPHA(I)*ALPHA(I)

1960 NEXT I

1970 FOR I=1 TO M

1980　 ALPHA(I)=ALPHA(I)/SQR(S)

1990 NEXT I

2000 '

2010 GOSUB *HYPERPLANE.DEF

2020 GOSUB *HYPERPLANE.1

2030 GOSUB *HYPERPLANE.2

2040 GOSUB *HYPERPLANE.3

2050 RETURN

2060 '

2070 *HYPERPLANE.DEF:

2080 M=M0

2090 KEY1$="123"

2100 KEY2$="456"

2110 MI=LEN(KEY1$)

2120 MJ=M0-MI

2130 RETURN

2140 '

2150 *HYPERPLANE.1:

2160 GOSUB *JACOBIAN

2170 GOSUB *HESSIAN

2180 '

2190 FOR I=1 TO M

2200　 FOR J=1 TO M

2210　　 X=HESSE(I,J)

2220　　 X1=INSTR(KEY1$,MID$(STR$(I),2))

2230　　 X2=INSTR(KEY2$,MID$(STR$(I),2))

2240　　 Y1=INSTR(KEY1$,MID$(STR$(J),2))

2250　　 Y2=INSTR(KEY2$,MID$(STR$(J),2))

2260　　 IF X1<>0　AND Y1<>0 THEN C1(X1,Y1)=X

2270　　 IF X1<>0　AND Y1=0　THEN C2(X1,Y2)=X

2280　　 IF X1=0　 AND Y1<>0 THEN C3(X2,Y1)=X

2290　　 IF X1=0　 AND Y1=0　THEN C4(X2,Y2)=X

2300　 NEXT J

2310 NEXT I

2320 RETURN

2330 '

2340 *HYPERPLANE.2:

2350 FOR I=1 TO M

2360　 X=ALPHA(I)

2370　 X1=INSTR(KEY1$,MID$(STR$(I),2))

2380　 X2=INSTR(KEY2$,MID$(STR$(I),2))

2390　 IF X1<>0 THEN ALPHA1(X1)=X

2400　 IF X2<>0 THEN ALPHA2(X2)=X

2410 NEXT I

2420 '

2430 ' [contour]

2440 FOR I=1 TO MJ

2450　 FOR J=1 TO MJ:AZ(I,J)=C4(I,J):NEXT J

2460 NEXT I

2470 NZ=MJ:GOSUB *INVERSE

2480 '

2490 FOR I=1 TO MJ

2500　 FOR J=1 TO MI

2510　　 S=0

2520　　 FOR K=1 TO MJ

2530　　　 S=S+BZ(I,K)*C3(K,J)

2540　　 NEXT K

2550　　 ALPHA4(I,J)=S

2560　 NEXT J

2570 NEXT I

2580 RETURN

2590 '

2600 *HYPERPLANE.3:

2610 RN=1:CN=MJ

2620 FOR J=1 TO CN:XX1(1,J)=ALPHA2(J):NEXT J

2630 GOSUB *PENROSE

2640 '

2650 FOR I=1 TO MJ

2660　 FOR J=1 TO MI

2670　　'ALPHA3(I,J)=YY1(I,1)*ALPHA1(J):' [arbitary]

2680　　'ALPHA3(I,J)=0　　　　　　　　 :' [section]

2690　　 ALPHA3(I,J)=ALPHA4(I,J)　　　 :' [contour]

2700　 NEXT J

2710 NEXT I

2720 '

2730 FOR I=1 TO MI

2740　 FOR J=1 TO MI

2750　　 S=0

2760　　 FOR K=1 TO MJ

2770　　　 S=S+C2(I,K)*ALPHA3(K,J)

2780　　 NEXT K

2790　　 C5(I,J)=S

2800　 NEXT J

2810 NEXT I

2820 '

2830 FOR I=1 TO MI

2840　 FOR J=1 TO MI

2850　　 S=0

2860　　 FOR K=1 TO MJ

2870　　　 S=S+ALPHA3(K,I)*C3(K,J)

2880　　 NEXT K

2890　　 C6(I,J)=S

2900　 NEXT J

2910 NEXT I

2920 '

2930 FOR I=1 TO MI

2940　 FOR J=1 TO MJ

2950　　 S=0

2960　　 FOR K=1 TO MJ

2970　　　 S=S+ALPHA3(K,I)*C4(K,J)

2980　　 NEXT K

2990　　 C7(I,J)=S

3000　 NEXT J

3010 NEXT I

3020 '

3030 FOR I=1 TO MI

3040　 FOR J=1 TO MI

3050　　 S=0

3060　　 FOR K=1 TO MJ

3070　　　 S=S+C7(I,K)*ALPHA3(K,J)

3080　　 NEXT K

3090　　 C8(I,J)=S

3100　 NEXT J

3110 NEXT I

3120 '

3130 FOR I=1 TO MI

3140　 FOR J=1 TO MI

3150　　'AZ(I,J)=C1(I,J)-C5(I,J)-C6(I,J)+C8(I,J)

3160　　 AZ(I,J)=C1(I,J)-C5(I,J)-C5(J,I)+C8(I,J)

3170　　 HESSE(I,J)=AZ(I,J)

3180　 NEXT J

3190 NEXT I

3200 NZ=MI:GOSUB *INVERSE

3210 '

3220 FOR I=1 TO MI

3230　 FOR J=1 TO MI:COV(I,J)=BZ(I,J):NEXT J

3240 NEXT I

3250 M0=M:M=MI

3260 RETURN

3270 '

3280 ' *** Moore-Penrose 逆行列 ( XX1==>YY1 ) ***

3290 '

3300 *PENROSE:

3310 GOSUB *EXACT2

3320 RETURN

3330 '

3340 ' ** XT*X **

3350 '

3360 *EXACT2:

3370 FOR I=1 TO CN

3380　 FOR J=I TO CN

3390　　 S=0

3400　　 FOR K=1 TO RN

3410　　　 S=S+XX1(K,I)*XX1(K,J)

3420　　 NEXT K

3430　　 AZ(I,J)=S:AZ(J,I)=S

3440　 NEXT J

3450 NEXT I

3460 '

3470 ' ** (XT*X)^(-1) **

3480 '

3490 NZ=CN:GOSUB *INVERSE

3500 '

3510 ' ** X*(XT*X)^(-1) **

3520 '

3530 FOR I=1 TO RN

3540　 FOR J=1 TO CN

3550　　 S=0

3560　　 FOR K=1 TO CN

3570　　　 S=S+XX1(I,K)*BZ(K,J)

3580　　 NEXT K

3590　　 H1(I,J)=S

3600　 NEXT J

3610 NEXT I

3620 '

3630 ' ** X*XT **

3640 '

3650 FOR I=1 TO RN

3660　 FOR J=I TO RN

3670　　 S=0

3680　　 FOR K=1 TO CN

3690　　　 S=S+XX1(I,K)*XX1(J,K)

3700　　 NEXT K

3710　　 AZ(I,J)=S:AZ(J,I)=S

3720　 NEXT J

3730 NEXT I

3740 '

3750 ' ** (X*XT)^(-1) **

3760 '

3770 NZ=RN:GOSUB *INVERSE

3780 '

3790 ' ** XT*(X*XT)^(-1) **

3800 '

3810 FOR I=1 TO CN

3820　 FOR J=1 TO RN

3830　　 S=0

3840　　 FOR K=1 TO RN

3850　　　 S=S+XX1(K,I)*BZ(K,J)

3860　　 NEXT K

3870　　 H2(I,J)=S

3880　 NEXT J

3890 NEXT I

3900 '

3910 ' ** XT*(X*XT)^(-1)*X*(XT*X)^(-1) **

3920 '

3930 FOR I=1 TO CN

3940　 FOR J=1 TO CN

3950　　 S=0

3960　　 FOR K=1 TO RN

3970　　　 S=S+H2(I,K)*H1(K,J)

3980　　 NEXT K

3990　　 H3(I,J)=S

4000　 NEXT J

4010 NEXT I

4020 '

4030 ' ** XT*(X*XT)^(-1)*X*(XT*X)^(-1)*XT **

4040 '

4050 FOR I=1 TO CN

4060　 FOR J=1 TO RN

4070　　 S=0

4080　　 FOR K=1 TO CN

4090　　　 S=S+H3(I,K)*XX1(J,K)

4100　　 NEXT K

4110　　 YY1(I,J)=S

4120　 NEXT J

4130 NEXT I

4140 RETURN

4150 '

4160 ' ***　超平面への正射影 ***

4170 '

4180 *HYPER.CONTOUR:

4190 GOSUB *HYPERPLANE.DEF

4200 GOSUB *HYPERPLANE.1

4210 '

4220 FOR I=1 TO MJ

4230　 FOR J=1 TO MJ:AZ(I,J)=C4(I,J):NEXT J

4240 NEXT I

4250 NZ=MJ:GOSUB *INVERSE

4260 '

4270 FOR I=1 TO MI

4280　 FOR J=1 TO MJ

4290　　 S=0

4300　　 FOR K=1 TO MJ

4310　　　 S=S+C2(I,K)*BZ(K,J)

4320　　 NEXT K

4330　　 C5(I,J)=S

4340　 NEXT J

4350 NEXT I

4360 '

4370 FOR I=1 TO MI

4380　 FOR J=1 TO MI

4390　　 S=0

4400　　 FOR K=1 TO MJ

4410　　　 S=S+C5(I,K)*C3(K,J)

4420　　 NEXT K

4430　　 C6(I,J)=S

4440　 NEXT J

4450 NEXT I

4460 '

4470 FOR I=1 TO MI

4480　 FOR J=1 TO MI

4490　　 AZ(I,J)=C1(I,J)-C6(I,J)

4500　　 HESSE(I,J)=AZ(I,J)

4510　 NEXT J

4520 NEXT I

4530 NZ=MI:GOSUB *INVERSE

4540 '

4550 FOR I=1 TO MI

4560　 FOR J=1 TO MI:COV(I,J)=BZ(I,J):NEXT J

4570 NEXT I

4580 M0=M:M=MI

4590 RETURN

4600 '

4610 ' ***　超平面による切断面 ***

4620 '

4630 *HYPER.SECTION:

4640 GOSUB *HYPERPLANE.DEF

4650 GOSUB *HYPERPLANE.1

4660 '

4670 FOR I=1 TO MI

4680　 FOR J=1 TO MI

4690　　 AZ(I,J)=C1(I,J)

4700　　 HESSE(I,J)=AZ(I,J)

4710　 NEXT J

4720 NEXT I

4730 NZ=MI:GOSUB *INVERSE

4740 '

4750 FOR I=1 TO MI

4760　 FOR J=1 TO MI:COV(I,J)=BZ(I,J):NEXT J

4770 NEXT I

4780 M0=M:M=MI

4790 RETURN

4800 '

4810 ' *** ヤコビアンの読み込み ***

4820 '

4830 *JACOBIAN:

4840 RESTORE *J

4850 FOR J=1 TO M

4860　 FOR I=1 TO N1

4870　　 READ JACOBI(J,I)

4880　 NEXT I

4890 NEXT J

4900 '

4910 RESTORE *V

4920 FOR J=1 TO M

4930　 READ V(J)

4940 NEXT J

4950 RETURN

4960 '

4970 *J:'[1]

4980 DATA　.2084, .0320, .0409,-.0413,-.0032,-.0086, .0033

4990 DATA　　　0,-.3443, .1829,-.1264,-.0153,-.0227,-.0493

5000 DATA -.3425, .0606,-.1273,-.1929, .0157, .0446,-.0244

5010 DATA　　　0, .8841, .4161, .1965, .9440,-.3250, .1699

5020 DATA　　　1,　　 0, .4550, .8003,-.3106,-.8159, .4470

5030 DATA　　　0,-.4673, .7873,-.5664,-.1115,-.4782,-.8782

5040 '

5050 *V:'[1]

5060 DATA -.0013

5070 DATA -.0020

5080 DATA -.0005

5090 DATA　.0007

5100 DATA -.0506

5110 DATA -.0119

5120 '

5130 '*J:'[20]

5140 DATA　.3063,-.1201, .1035, .1037,-.0128,-.0019, .0373

5150 DATA　　　0,-.1692, .2014,-.2040,-.0076,-.0405,-.0248

5160 DATA -.1289, .1619,-.1761,-.1737, .0217, .0301,-.0311

5170 DATA　　　0, .8031,-.4195, .3965, .8705, .3581, .0393

5180 DATA　　　1,　　 0, .7101, .7037,-.0844,-.5671, .8029

5190 DATA　　　0, .5958, .5655,-.5896, .4848,-.7418,-.5948

5200 '

5210 '*V:'[20]

5220 DATA　-.6581

5230 DATA　 .0862

5240 DATA　 .1687

5250 DATA　 .5254

5260 DATA -1.1528

5270 DATA　 .7858

5280 '

5290 ' *** ヘシアンと共分散行列 ***

5300 '

5310 *HESSIAN:

5320 FOR M1=1 TO M

5330　 FOR M2=M1 TO M

5340　　 S=0

5350　　 FOR I=1 TO N1

5360　　　 S=S+JACOBI(M1,I)*JACOBI(M2,I)

5370　　 NEXT I

5380　　 AZ(M1,M2)=S:AZ(M2,M1)=S

5390　　 COV(M1,M2)=S:COV(M2,M1)=S

5400　 NEXT M2

5410 NEXT M1

5420 NZ=M

5430 '

5440 GOSUB *INVERSE

5450 '

5460 FOR I=1 TO M

5470　 FOR J=1 TO M:HESSE(I,J)=BZ(I,J):NEXT J

5480 NEXT I

5490 RETURN

5500 '

5510 ' ** 逆行列 ( AZ==>BZ ) **

5520 '

5530 *INVERSE:' [GAUSS-JORDAN]

5540 FOR I=1 TO NZ

5550　 FOR J=1 TO NZ:BZ(I,J)=0:NEXT J

5560　 BZ(I,I)=1

5570 NEXT I

5580 FOR I=1 TO NZ

5590　 IZ=I

5600　 IF AZ(I,I)=0 THEN 5610 ELSE 5670

5610　 IZ=IZ+1:IF IZ>NZ THEN RZ=1:RETURN

5620　 IF AZ(IZ,I)=0 THEN 5610

5630　 FOR J=1 TO NZ

5640　　 SWAP AZ(I,J),AZ(IZ,J)

5650　　 SWAP BZ(I,J),BZ(IZ,J)

5660　 NEXT J

5670　 AII=AZ(I,I)

5680　 FOR J=1 TO NZ

5690　　 AZ(I,J)=AZ(I,J)/AII

5700　　 BZ(I,J)=BZ(I,J)/AII

5710　 NEXT J

5720　 FOR K=1 TO NZ

5730　　 AKI=AZ(K,I):IF K=I THEN 5780

5740　　 FOR J=1 TO NZ

5750　　　 AZ(K,J)=AZ(K,J)-AKI*AZ(I,J)

5760　　　 BZ(K,J)=BZ(K,J)-AKI*BZ(I,J)

5770　　 NEXT J

5780　 NEXT K

5790 NEXT I

5800 RZ=0

5810 RETURN

5820 '

5830 ' *** 楕円体 ***

5840 '

5850 *ELLIPSOID:

5860 IF M<2 THEN RETURN

5870 '

5880 FOR I=1 TO M

5890　 FOR J=1 TO M:AZ(I,J)=COV(I,J):NEXT J

5900 NEXT I

5910 NZ=M:GOSUB *JACOBI

5920 '

5930 FOR I=1 TO M-1

5940　 FOR J=I+1 TO M

5950　　 IB=VAL(MID$(KEY1$,I,1))

5960　　 JB=VAL(MID$(KEY1$,J,1))

5970　　 UO=VO(IB):UL=VL(IB)

5980　　 VO=VO(JB):VL=VL(JB)

5990　　 DU=VD(IB):SU=VS(IB)

6000　　 DV=VD(JB):SV=VS(JB)

6010 '

6020　　 GOSUB *DRAW.AXIS

6030　　 GOSUB *DRAW.ELLIPSE:'IF I=1 AND J=2 THEN FL$="b:e1":GOSUB *GSAVE

6040　　 GOSUB *TEMP

6050　 NEXT J

6060 NEXT I

6070 RETURN

6080 '

6090 ' ** 固有値 ( AZ==>BZ ) **

6100 '

6110 *JACOBI:' [ EIGEN-VECTOR: BZ(J,I) ]

6120 EPSL=.00001

6130 FOR I=1 TO NZ

6140　 FOR J=1 TO NZ:BZ(I,J)=0:NEXT J

6150　 BZ(I,I)=1

6160 NEXT I

6170 '

6180 CMAX=0

6190 FOR I=1 TO NZ

6200　 FOR J=1 TO NZ

6210　　 IF CMAX<ABS(AZ(I,J)) THEN CMAX=ABS(AZ(I,J))

6220　 NEXT J

6230 NEXT I

6240 EPS=CMAX*EPSL

6250 '

6260 LMAX=100

6270 LOOP=0

6280 SW=0

6290 WHILE LOOP<LMAX AND SW=0

6300 '

6310 FOR L=1 TO NZ-1

6320 FOR I=1 TO NZ-1

6330　 J=I+L:IF J>NZ THEN 6520

6340　 IF ABS(AZ(I,J))<EPS THEN 6510

6350　 DIF=AZ(I,I)-AZ(J,J)

6360　 IF DIF=0 THEN T=PI/4 ELSE T=1/2*ATN(2*AZ(I,J)/DIF)

6370　 C=COS(T):S=SIN(T)

6380　 FOR K=1 TO NZ

6390　　 A1=AZ(I,K):A2=AZ(J,K)

6400　　 AZ(I,K)= C*A1+S*A2

6410　　 AZ(J,K)=-S*A1+C*A2

6420　 NEXT K

6430　 FOR K=1 TO NZ

6440　　 A1=AZ(K,I):A2=AZ(K,J)

6450　　 AZ(K,I)= C*A1+S*A2

6460　　 AZ(K,J)=-S*A1+C*A2

6470　　 B1=BZ(K,I):B2=BZ(K,J)

6480　　 BZ(K,I)= C*B1+S*B2

6490　　 BZ(K,J)=-S*B1+C*B2

6500　 NEXT K

6510 NEXT I

6520 NEXT L

6530 '

6540 SW=1

6550 FOR I=1 TO NZ-1

6560　 FOR J=I+1 TO NZ

6570　　 IF ABS(AZ(I,J))>EPS THEN SW=0

6580　 NEXT J

6590 NEXT I

6600 LOOP=LOOP+1

6610 '

6620 WEND

6630 RZ=0

6640 RETURN

6650 '

6660 ' *** 正射影 ***

6670 '

6680 *DRAW.ELLIPSE:

6690 'DRAW.COLOR=4

6700 WUO=UO-UL:WUL=UO+UL

6710 WVO=VO-VL:WVL=VO+VL

6720 WINDOW(WUO,-WVL)-(WUL,-WVO)

6730 VIEW(SX1,SY1)-(SX2,SY2)

6740 '

6750 'GOSUB *POLAR

6760 'GOSUB *EQUATOR

6770 'GOSUB *SECTION

6780 GOSUB *CONTOUR

6790 'GOSUB *HYPERCUBE

6800 'GOSUB *DUALCUBE

6810 'GOSUB *SIMPLEX

6820 '

6830 WINDOW(0,0)-(639,399)

6840 VIEW(0,0)-(639,399)

6850 RETURN

6860 '

6870 ' *** 輪郭 ***

6880 '

6890 *CONTOUR:

6900 GOTO 7120

6910 '

6920 S1=COV(I,I)

6930 S2=COV(J,J)

6940 S3=COV(I,J)

6950 LAMBDA1=(S1+S2+SQR((S1-S2)*(S1-S2)+4*S3*S3))/2

6960 LAMBDA2=(S1+S2-SQR((S1-S2)*(S1-S2)+4*S3*S3))/2

6970 'LAMBDA2=(S1*S2-S3*S3)/LAMBDA1

6980 T=1:' [non-stochastic]

6990 AAA=SQR(T*LAMBDA1)

7000 BBB=SQR(T*LAMBDA2)

7010 IF S3=0 THEN CCC=PI/2　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 ELSE B3=(S2-S1+SQR((S1-S2)*(S1-S2)+4*S3*S3))/2/S3:CCC=ATN(B3)

7020 'IF S3=0 THEN CCC=PI/2 ELSE B3=(LAMBDA1-S2)/S3:CCC=ATN(B3)

7030 '

7040 G=0

7050 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

7060　 GX=AAA*COS(ANGLE)*COS(CCC)-BBB*SIN(ANGLE)*SIN(CCC)+VO(IB)

7070　 GY=AAA*COS(ANGLE)*SIN(CCC)+BBB*SIN(ANGLE)*COS(CCC)+VO(JB)

7080　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

7090　 LINE-(GX,-GY),DRAW.COLOR

7100 NEXT ANGLE

7110 RETURN

7120 '

7130 DIF=COV(I,I)-COV(J,J)

7140 IF DIF=0 THEN TH=PI/4 ELSE TH=1/2*ATN(2*COV(I,J)/DIF)

7150 C=COS(TH):S=SIN(TH)

7160 '

7170 LAMBDA1=COV(I,I)*C*C+2*COV(I,J)*C*S+COV(J,J)*S*S

7180 LAMBDA2=COV(I,I)*S*S-2*COV(I,J)*C*S+COV(J,J)*C*C

7190 T=1:' [non-stochastic]

7200 AAA=SQR(T*LAMBDA1)

7210 BBB=SQR(T*LAMBDA2)

7220 '

7230 G=0

7240 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

7250　 GX=AAA*COS(ANGLE)*C-BBB*SIN(ANGLE)*S+VO(IB)

7260　 GY=AAA*COS(ANGLE)*S+BBB*SIN(ANGLE)*C+VO(JB)

7270　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

7280　 LINE-(GX,-GY),DRAW.COLOR

7290 NEXT ANGLE

7300 RETURN

7310 '

7320 ' *** 極線 ***

7330 '

7340 *POLAR:

7350 FOR K=1 TO M

7360　 T=1:' [non-stochastic]

7370　 RRR=SQR(T*AZ(K,K))

7380　 GX=BZ(I,K)*RRR

7390　 GY=BZ(J,K)*RRR

7400　 GX1= GX+VO(IB):GY1= GY+VO(JB)

7410　 GX2=-GX+VO(IB):GY2=-GY+VO(JB)

7420　 LINE(GX1,-GY1)-(GX2,-GY2),7

7430　 IF K=I OR K=J THEN LINE(GX1,-GY1)-(GX2,-GY2),CLR

7440 NEXT K

7450 RETURN

7460 '

7470 ' *** 赤道 ***

7480 '

7490 *EQUATOR:

7500 T=1:' [non-stochastic]

7510 AAA=SQR(T*AZ(I,I))

7520 BBB=SQR(T*AZ(J,J))

7530 '

7540 G=0

7550 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

7560　 GX=AAA*COS(ANGLE)*BZ(I,I)+BBB*SIN(ANGLE)*BZ(I,J)+VO(IB)

7570　 GY=AAA*COS(ANGLE)*BZ(J,I)+BBB*SIN(ANGLE)*BZ(J,J)+VO(JB)

7580　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

7590　 LINE-(GX,-GY),DRAW.COLOR

7600 NEXT ANGLE

7610 RETURN

7620 '

7630 ' *** 断面 ***

7640 '

7650 *SECTION:

7660 DIF=HESSE(I,I)-HESSE(J,J)

7670 IF DIF=0 THEN TH=PI/4 ELSE TH=1/2*ATN(2*HESSE(I,J)/DIF)

7680 C=COS(TH):S=SIN(TH)

7690 '

7700 LAMBDA1=HESSE(I,I)*C*C+2*HESSE(I,J)*C*S+HESSE(J,J)*S*S

7710 LAMBDA2=HESSE(I,I)*S*S-2*HESSE(I,J)*C*S+HESSE(J,J)*C*C

7720 T=1:' [non-stochastic]

7730 AAA=SQR(T/LAMBDA1)

7740 BBB=SQR(T/LAMBDA2)

7750 '

7760 G=0

7770 FOR ANGLE=0 TO 2*PI+.01 STEP PI/180

7780　 GX=AAA*COS(ANGLE)*C-BBB*SIN(ANGLE)*S+VO(IB)

7790　 GY=AAA*COS(ANGLE)*S+BBB*SIN(ANGLE)*C+VO(JB)

7800　 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

7810　 LINE-(GX,-GY),DRAW.COLOR

7820 NEXT ANGLE

7830 RETURN

7840 '

7850 ' *** 超立方体 ***

7860 '

7870 *HYPERCUBE:

7880 T=1:' [non-stochastic]

7890 FOR IS=1 TO M:EDGE(IS)=SQR(T*AZ(IS,IS)):NEXT IS

7900 '

7910 FOR IS=0 TO 2^M-1

7920　 FOR JS=1 TO M:S(JS)=-1:NEXT JS

7930　 FOR JS=1 TO M

7940　　 IF (IS AND 2^(JS-1))=0 THEN S(JS)=1

7950　 NEXT JS

7960　 FOR JS=1 TO M:S0(JS)=S(JS):NEXT JS

7970　'GOSUB *VERTEX:GXS=GX:GYS=GY

7980　 '

7990　 COUNT=0

8000　 FOR JS=1 TO M

8010　　 COUNT=COUNT+1

8020　　 FOR KS=1 TO M:S(KS)=S0(KS):NEXT KS

8030　　 S(JS)=-S(JS)

8040　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

8050　　'GOSUB *VERTEX:GXE=GX:GYE=GY

8060　　'LINE(GXS,-GYS)-(GXE,-GYE),7

8070　 NEXT JS

8080　 GOSUB *POLYGON

8090 NEXT IS

8100 RETURN

8110 '

8120 ' *** 頂点 ***

8130 '

8140 *VERTEX:

8150 'GOSUB *INTERNAL

8160 GX=0:GY=0

8170 FOR GS=1 TO M

8180　 GX=GX+BZ(I,GS)*EDGE(GS)*S(GS)

8190　 GY=GY+BZ(J,GS)*EDGE(GS)*S(GS)

8200 NEXT GS

8210 GX=GX+VO(IB):GY=GY+VO(JB)

8220 RETURN

8230 '

8240 ' *** 内接化 ***

8250 '

8260 *INTERNAL:

8270 X=(1-SQR(1+M))/M

8280 G=SQR(1+1/M)

8290 FOR GS=1 TO M

8300　 S(GS)=(S(GS)-(X+1)/(M+1))*G

8310 NEXT GS

8320 RETURN

8330 '

8340 ' *** 輪郭抽出 ***

8350 '

8360 '*POLYGON:

8370 EPS=.05

8380 FOR KS=1 TO M

8390　 FOR LS=1 TO M:S(LS)=S0(LS):NEXT LS

8400　 IF KS=JS THEN 8420

8410　 S(KS)=S(KS)*(1-EPS)

8420　 GOSUB *VERTEX:GXS(KS)=GX:GYS(KS)=GY

8430　 '

8440　 FOR LS=1 TO M:S(LS)=S0(LS):NEXT LS

8450　 S(JS)=-S(JS)

8460　 IF KS=JS THEN 8480

8470　 S(KS)=S(KS)*(1-EPS)

8480　 GOSUB *VERTEX:GXE(KS)=GX:GYE(KS)=GY

8490 NEXT KS

8500 '

8510 X=GXS(JS):Y=GYS(JS)

8520 FOR KS=1 TO M

8530　 GZ(KS)=(GYE(KS)-GYS(KS))*(X-GXS(KS))-(GXE(KS)-GXS(KS))*(Y-GYS(KS))

8540 NEXT KS

8550 '

8560 ZMAX=GZ(JS):ZMIN=GZ(JS)

8570 FOR KS=1 TO M

8580　 IF KS=JS THEN 8610

8590　 IF ZMAX<GZ(KS) THEN ZMAX=GZ(KS)

8600　 IF ZMIN>GZ(KS) THEN ZMIN=GZ(KS)

8610 NEXT KS

8620 '

8630 SW=0

8640 IF (ZMAX-GZ(JS))*(ZMIN-GZ(JS))=0 THEN SW=1

8650 GXS=GXS(JS):GYS=GYS(JS)

8660 GXE=GXE(JS):GYE=GYE(JS)

8670 IF SW=1 THEN LINE(GXS,-GYS)-(GXE,-GYE),7

8680 RETURN

8690 '

8700 ' *** 輪郭抽出 ***

8710 '

8720 *POLYGON:

8730 EPS=.05

8740 FOR JS=1 TO COUNT

8750　 FOR LS=1 TO M:S1(LS)=SS(JS,LS):NEXT LS

8760　 FOR KS=1 TO COUNT

8770　　 FOR LS=1 TO M:S2(LS)=SS(KS,LS):NEXT LS

8780　　 '

8790　　 FOR LS=1 TO M:S(LS)=S0(LS):NEXT LS

8800　　 IF KS=JS THEN 8820

8810　　 FOR LS=1 TO M:S(LS)=S0(LS)+EPS*(S2(LS)-S0(LS)):NEXT LS

8820　　 GOSUB *VERTEX:GXS(KS)=GX:GYS(KS)=GY

8830　　 '

8840　　 FOR LS=1 TO M:S(LS)=S1(LS):NEXT LS

8850　　 IF KS=JS THEN 8870

8860　　 FOR LS=1 TO M:S(LS)=S1(LS)+EPS*(S2(LS)-S0(LS)):NEXT LS

8870　　 GOSUB *VERTEX:GXE(KS)=GX:GYE(KS)=GY

8880　 NEXT KS

8890　 GOSUB *POLYGON.LINE

8900 NEXT JS

8910 RETURN

8920 '

8930 *POLYGON.LINE:

8940 X=GXS(JS):Y=GYS(JS)

8950 FOR KS=1 TO COUNT

8960　 GZ(KS)=(GYE(KS)-GYS(KS))*(X-GXS(KS))-(GXE(KS)-GXS(KS))*(Y-GYS(KS))

8970 NEXT KS

8980 '

8990 ZMAX=GZ(JS):ZMIN=GZ(JS)

9000 FOR KS=1 TO COUNT

9010　 IF KS=JS THEN 9040

9020　 IF ZMAX<GZ(KS) THEN ZMAX=GZ(KS)

9030　 IF ZMIN>GZ(KS) THEN ZMIN=GZ(KS)

9040 NEXT KS

9050 '

9060 SW=0

9070 IF (ZMAX-GZ(JS))*(ZMIN-GZ(JS))=0 THEN SW=1

9080 GXS=GXS(JS):GYS=GYS(JS)

9090 GXE=GXE(JS):GYE=GYE(JS)

9100 IF SW=1 THEN LINE(GXS,-GYS)-(GXE,-GYE),7

9110 RETURN

9120 '

9130 ' *** 双対立方体 ***

9140 '

9150 *DUALCUBE:

9160 T=1:' [non-stochastic]

9170 FOR IS=1 TO M:EDGE(IS)=SQR(T*AZ(IS,IS)):NEXT IS

9180 '

9190 FOR IS=1 TO M

9200　 FOR JS=1 TO M:S(JS)=0:NEXT JS

9210　 S(IS)=1

9220　 FOR JS=1 TO M:S0(JS)=S(JS):NEXT JS

9230　'GOSUB *VERTEX:GXS=GX:GYS=GY

9240　 '

9250　 COUNT=0

9260　 FOR JS=1 TO M

9270　　 IF JS=IS THEN 9410

9280　　 COUNT=COUNT+1

9290　　 FOR KS=1 TO M:S(KS)=0:NEXT KS

9300　　 S(JS)=1

9310　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9320　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9330　　'LINE(GXS,-GYS)-(GXE,-GYE),7

9340　　 '

9350　　 COUNT=COUNT+1

9360　　 FOR KS=1 TO M:S(KS)=0:NEXT KS

9370　　 S(JS)=-1

9380　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9390　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9400　　'LINE(GXS,-GYS)-(GXE,-GYE),7

9410　 NEXT JS

9420　 GOSUB *POLYGON

9430 NEXT IS

9440 '

9450 FOR IS=1 TO M

9460　 FOR JS=1 TO M:S(JS)=0:NEXT JS

9470　 S(IS)=-1

9480　 FOR JS=1 TO M:S0(JS)=S(JS):NEXT JS

9490　'GOSUB *VERTEX:GXS=GX:GYS=GY

9500　 '

9510　 COUNT=0

9520　 FOR JS=1 TO M

9530　　 IF JS=IS THEN 9670

9540　　 COUNT=COUNT+1

9550　　 FOR KS=1 TO M:S(KS)=0:NEXT KS

9560　　 S(JS)=1

9570　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9580　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9590　　'LINE(GXS,-GYS)-(GXE,-GYE),7

9600　　 '

9610　　 COUNT=COUNT+1

9620　　 FOR KS=1 TO M:S(KS)=0:NEXT KS

9630　　 S(JS)=-1

9640　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9650　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9660　　'LINE(GXS,-GYS)-(GXE,-GYE),7

9670　 NEXT JS

9680　 GOSUB *POLYGON

9690 NEXT IS

9700 RETURN

9710 '

9720 ' *** 正単体 ***

9730 '

9740 *SIMPLEX:

9750 T=1:' [non-stochastic]

9760 FOR IS=1 TO M:EDGE(IS)=SQR(T*AZ(IS,IS)):NEXT IS

9770 '

9780 FOR IS=1 TO M

9790　 FOR JS=1 TO M:S(JS)=0:NEXT JS

9800　 S(IS)=1

9810　 FOR JS=1 TO M:S0(JS)=S(JS):NEXT JS

9820　'GOSUB *VERTEX:GXS=GX:GYS=GY

9830　 '

9840　 COUNT=0

9850　 FOR JS=1 TO M

9860　　 IF JS=IS THEN 9930

9870　　 COUNT=COUNT+1

9880　　 FOR KS=1 TO M:S(KS)=0:NEXT KS

9890　　 S(JS)=1

9900　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9910　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9920　　'LINE(GXS,-GYS)-(GXE,-GYE),7

9930　 NEXT JS

9940　 '

9950　　 COUNT=COUNT+1

9960　　 FOR KS=1 TO M:S(KS)=(1-SQR(1+M))/M:NEXT KS

9970　　 FOR KS=1 TO M:SS(COUNT,KS)=S(KS):NEXT KS

9980　　'GOSUB *VERTEX:GXE=GX:GYE=GY

9990　　'LINE(GXS,-GYS)-(GXE,-GYE),7

10000　 GOSUB *POLYGON

10010 NEXT IS

10020 RETURN

10030 '

10040 ' *** 座標軸 ***

10050 '

10060 *DRAW.AXIS:

10070 CLS 3

10080 SX1=59:SY1=360*(1-RY):SX2=580*RX+59:SY2=360

10090 LINE(SX1,SY1)-(SX2,SY2),GRP.COLOR,B

10100 RR=3:GOSUB *X.SCALE2:GOSUB *X.REF2

10110 RR=3:GOSUB *Y.SCALE2:GOSUB *Y.REF2

10120 RR=3:GOSUB *DOT.PLOT

10130 GOSUB *TITLE.BACK2

10140 RETURN

10150 '

10160 ' *** 等分目盛り (X) ***

10170 '

10180 *X.SCALE2:

10190 CU=(SX2-SX1)/DU/2

10200 IF KIND=0 THEN KARA=-RR:MADE=0

10210 IF KIND=1 THEN KARA=0　:MADE=RR

10220 IF KIND=2 THEN KARA=-RR:MADE=RR

10230 FOR K=0 TO DU*2

10240　 LINE(CU*K+SX1,SY2+KARA)-(CU*K+SX1,SY2+MADE),GRP.COLOR

10250　 LINE(CU*K+SX1,SY1-KARA)-(CU*K+SX1,SY1-MADE),GRP.COLOR

10260 NEXT K

10270 'GOSUB *X.REF2

10280 RETURN

10290 '

10300 ' *** 参照値の記入 (X) ***

10310 '

10320 *X.REF2:

10330 REF.COLOR=6

10340 UW=UL-UO

10350 FOR K=0 TO DU*2 STEP SU

10360　 F=UW*(K/DU-1):F$=STR$(F)

10370　 IF F>=0 THEN F$=MID$(F$,2)

10380　 FL=LEN(F$)

10390　 FOR L=1 TO FL

10400　　 GX=CU*K+59-FL*8/2+(L-1)*8:GY=23*16

10410　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

10420　 NEXT L

10430 NEXT K

10440 RETURN

10450 '

10460 ' *** 等分目盛り (Y) ***

10470 '

10480 *Y.SCALE2:

10490 CV=(SY2-SY1)/DV/2

10500 IF KIND=0 THEN KARA=0　:MADE=RR

10510 IF KIND=1 THEN KARA=-RR:MADE=0

10520 IF KIND=2 THEN KARA=-RR:MADE=RR

10530 FOR K=0 TO DV*2

10540　 LINE(SX1+KARA,SY2-CV*K)-(SX1+MADE,SY2-CV*K),GRP.COLOR

10550　 LINE(SX2-KARA,SY2-CV*K)-(SX2-MADE,SY2-CV*K),GRP.COLOR

10560 NEXT K

10570 'GOSUB *Y.REF2

10580 RETURN

10590 '

10600 ' *** 参照値の記入 (Y) ***

10610 '

10620 *Y.REF2:

10630 REF.COLOR=6

10640 VW=VL-VO

10650 FOR K=0 TO DV*2 STEP SV

10660　 F=VW*(K/DV-1):F$=STR$(F)

10670　 IF F>=0 THEN F$=MID$(F$,2)

10680　 FL=LEN(F$)

10690　 FOR L=1 TO FL

10700　　 GX=48-FL*8+(L-1)*8:GY=360-CV*K-8

10710　　 PUT(GX,GY),KANJI(ASC(MID$(F$,L,1))),PSET,REF.COLOR,0

10720　 NEXT L

10730 NEXT K

10740 RETURN

10750 '

10760 ' *** プロット ***

10770 '

10780 *DOT.PLOT:

10790 WUO=UO-UL:WUL=UO+UL

10800 WVO=VO-VL:WVL=VO+VL

10810 WX=WUL-WUO

10820 BX=(WUL*SX1-WUO*SX2)/WX

10830 AX=(SX2-SX1)/WX

10840 WY=WVL-WVO

10850 BY=(WVL*SY2-WVO*SY1)/WY

10860 AY=(SY1-SY2)/WY

10870 '

10880 XX=AX*V(IB)+BX

10890 YY=AY*V(JB)+BY

10900 CIRCLE(XX,YY),RR,1

10910 PAINT(XX,YY),CLR,1

10920 CIRCLE(XX,YY),RR,CLR

10930 RETURN

10940 '

10950 ' *** タイトル ***

10960 '

10970 *TITLE.BACK2:

10980 LX=INT((580*RX+59)/8)

10990 LY=INT(360*(1-RY)/16)

11000 LZ=INT((290*RX+59)/8)

11010 L1$="ellipsoid projection of "

11020 L2$="v"+MID$(STR$(IB),2)+"-"

11030 L3$="v"+MID$(STR$(JB),2)+" plane"

11040 L$=L1$+L2$+L3$

11050 LOCATE LZ-LEN(L$)\2,LY-3:PRINT L$

11060 L$=UNIT$(JB)

11070 IF LEN(L$)=0 THEN 11090

11080 LOCATE 7,LY-1:PRINT "("+L$+")"

11090 L$=UNIT$(IB)

11100 IF LEN(L$)=0 THEN 11120

11110 LOCATE LX+3,21:PRINT "("+L$+")"

11120 '

11130 LOCATE LX+3,LY:PRINT "v components: "

11140 FOR K=1 TO M0

11150　 V$="v"+MID$(STR$(K),2)+"="

11160　 LOCATE LX+3,LY+K+1:PRINT V$;:PRINT V(K)

11170 NEXT K

11180 RETURN

11190 '

11200 *TEMP:

11210 LOCATE 0,24:PRINT "何かキーを押してください";:WHILE INKEY$="":WEND

11220 RETURN

11230 '

11240 ' *** ワイヤーフレーム楕円体 ***

11250 '

11260 *ELLIPSOID.WF:

11270 IF M<3 THEN RETURN

11280 '

11290 IB=VAL(MID$(KEY1$,1,1))

11300 JB=VAL(MID$(KEY1$,2,1))

11310 KB=VAL(MID$(KEY1$,3,1))

11320 IJK=3:GOSUB *ELLIPSOID.WF.SUB

11330 '

11340 IB=VAL(MID$(KEY1$,1,1))

11350 JB=VAL(MID$(KEY1$,3,1))

11360 KB=VAL(MID$(KEY1$,2,1))

11370 IJK=3:GOSUB *ELLIPSOID.WF.SUB

11380 '

11390 IB=VAL(MID$(KEY1$,2,1))

11400 JB=VAL(MID$(KEY1$,3,1))

11410 KB=VAL(MID$(KEY1$,1,1))

11420 IJK=3:GOSUB *ELLIPSOID.WF.SUB

11430 RETURN

11440 '

11450 *ELLIPSOID.WF.SUB:

11460 UO=VO(IB):UL=VL(IB)

11470 VO=VO(JB):VL=VL(JB)

11480 DU=VD(IB):SU=VS(IB)

11490 DV=VD(JB):SV=VS(JB)

11500 '

11510 GOSUB *DRAW.AXIS.WF

11520 GOSUB *DRAW.ELLIPSE.WF:'IF IB=1 AND JB=2 THEN FL$="b:e2":GOSUB *GSAVE

11530 GOSUB *TEMP

11540 RETURN

11550 '

11560 ' *** 正射影 ***

11570 '

11580 *DRAW.ELLIPSE.WF:

11590 'DRAW.COLOR=4

11600 WUO=UO-UL:WUL=UO+UL

11610 WVO=VO-VL:WVL=VO+VL

11620 WINDOW(WUO,-WVL)-(WUL,-WVO)

11630 VIEW(SX1,SY1)-(SX2,SY2)

11640 '

11650 IF IJK<>0 THEN GOSUB *EIGEN.WF:GOSUB *COMMON.WF

11660 IF IJK=0 AND STAGE$="1" THEN I=IB:J=JB:GOSUB *CONTOUR

11670 IF IJK=0 AND STAGE$="2" THEN I=IB:J=JB:GOSUB *SECTION

11680 'GOSUB *COMMON.WF

11690 'I=IB:J=JB:GOSUB *CONTOUR

11700 'I=IB:J=JB:GOSUB *SECTION

11710 '

11720 WINDOW(0,0)-(639,399)

11730 VIEW(0,0)-(639,399)

11740 RETURN

11750 '

11760 *EIGEN.WF:

11770 IF STAGE$="2" THEN 11960

11780 '

11790 AZ(1,1)=COV(IB,IB)

11800 AZ(1,2)=COV(IB,JB)

11810 AZ(1,3)=COV(IB,KB)

11820 AZ(2,1)=AZ(1,2)

11830 AZ(2,2)=COV(JB,JB)

11840 AZ(2,3)=COV(JB,KB)

11850 AZ(3,1)=AZ(1,3)

11860 AZ(3,2)=AZ(2,3)

11870 AZ(3,3)=COV(KB,KB)

11880 NZ=3:GOSUB *JACOBI

11890 '

11900 T=1:' [non-stochastic]

11910 AAA=SQR(T*AZ(1,1))

11920 BBB=SQR(T*AZ(2,2))

11930 CCC=SQR(T*AZ(3,3))

11940 KS=PI/12

11950 GOTO 12130

11960 '

11970 AZ(1,1)=HESSE(IB,IB)

11980 AZ(1,2)=HESSE(IB,JB)

11990 AZ(1,3)=HESSE(IB,KB)

12000 AZ(2,1)=AZ(1,2)

12010 AZ(2,2)=HESSE(JB,JB)

12020 AZ(2,3)=HESSE(JB,KB)

12030 AZ(3,1)=AZ(1,3)

12040 AZ(3,2)=AZ(2,3)

12050 AZ(3,3)=HESSE(KB,KB)

12060 NZ=3:GOSUB *JACOBI

12070 '

12080 T=1:' [non-stochastic]

12090 AAA=SQR(T/AZ(1,1))

12100 BBB=SQR(T/AZ(2,2))

12110 CCC=SQR(T/AZ(3,3))

12120 KS=PI/12

12130 RETURN

12140 '

12150 ' *** 経線・緯線 ***

12160 '

12170 *COMMON.WF:

12180 FOR V=PI/2-KS TO -PI/2+KS-.01 STEP -KS

12190　 G=0

12200　 FOR U=0 TO 2*PI+.01 STEP KS

12210　　 GOSUB *WIRE.FRAME

12220　 NEXT U

12230 NEXT V

12240 '

12250 FOR U=0 TO 2*PI+.01 STEP KS

12260　 G=0

12270　'FOR V=PI/2-KS TO -PI/2+KS-.01 STEP -KS

12280　 FOR V=PI/2 TO -PI/2-.01 STEP -KS

12290　　 GOSUB *WIRE.FRAME

12300　 NEXT V

12310 NEXT U

12320 RETURN

12330 '

12340 *WIRE.FRAME:

12350 ON IJK GOTO 12480,12420,12360

12360 '

12370 '[k-pole]

12380 GX=AAA*COS(U)*COS(V)*BZ(1,1)+BBB*SIN(U)*COS(V)*BZ(1,2)+CCC*SIN(V)*BZ(1,3)

12390 GY=AAA*COS(U)*COS(V)*BZ(2,1)+BBB*SIN(U)*COS(V)*BZ(2,2)+CCC*SIN(V)*BZ(2,3)

12400 GZ=AAA*COS(U)*COS(V)*BZ(3,1)+BBB*SIN(U)*COS(V)*BZ(3,2)+CCC*SIN(V)*BZ(3,3)

12410 GOTO 12530

12420 '

12430 '[j-pole]

12440 GX=AAA*SIN(U)*COS(V)*BZ(1,1)+BBB*SIN(V)*BZ(1,2)+CCC*COS(U)*COS(V)*BZ(1,3)

12450 GY=AAA*SIN(U)*COS(V)*BZ(2,1)+BBB*SIN(V)*BZ(2,2)+CCC*COS(U)*COS(V)*BZ(2,3)

12460 GZ=AAA*SIN(U)*COS(V)*BZ(3,1)+BBB*SIN(V)*BZ(3,2)+CCC*COS(U)*COS(V)*BZ(3,3)

12470 GOTO 12530

12480 '

12490 '[i-pole]

12500 GX=AAA*SIN(V)*BZ(1,1)+BBB*COS(U)*COS(V)*BZ(1,2)+CCC*SIN(U)*COS(V)*BZ(1,3)

12510 GY=AAA*SIN(V)*BZ(2,1)+BBB*COS(U)*COS(V)*BZ(2,2)+CCC*SIN(U)*COS(V)*BZ(2,3)

12520 GZ=AAA*SIN(V)*BZ(3,1)+BBB*COS(U)*COS(V)*BZ(3,2)+CCC*SIN(U)*COS(V)*BZ(3,3)

12530 '

12540 GX=GX+VO(IB)

12550 GY=GY+VO(JB)

12560 GZ=GZ+VO(KB)

12570 '

12580 ON VHP+1 GOSUB *VERTICAL,*VERTICAL,*HORIZONTAL,*PROFILE

12590 RETURN

12600 '

12610 *VERTICAL:

12620 IF G=0 THEN PSET(GX,-GY),DRAW.COLOR:G=1

12630 GOSUB *NORMAL.VECTOR

12640 IF DZ<=0 THEN LINE-(GX,-GY),DRAW.COLOR

12650 IF DZ> 0 THEN LINE-(GX,-GY),DRAW.COLOR,,&H1111

12660 RETURN

12670 '

12680 *HORIZONTAL:

12690 IF G=0 THEN PSET(GX,-GZ),DRAW.COLOR:G=1

12700 GOSUB *NORMAL.VECTOR:DY=-DY

12710 IF DY<=0 THEN LINE-(GX,-GZ),DRAW.COLOR

12720 IF DY> 0 THEN LINE-(GX,-GZ),DRAW.COLOR,,&H1111

12730 RETURN

12740 '

12750 *PROFILE:

12760 IF G=0 THEN PSET(GZ,-GY),DRAW.COLOR:G=1

12770 GOSUB *NORMAL.VECTOR:DX=-DX

12780 IF DX<=0 THEN LINE-(GZ,-GY),DRAW.COLOR

12790 IF DX> 0 THEN LINE-(GZ,-GY),DRAW.COLOR,,&H1111

12800 RETURN

12810 '

12820 ' *** 法線ベクトル ***

12830 '

12840 *NORMAL.VECTOR:

12850 ON IJK GOTO 13040,12950,12860

12860 '

12870 '[k-pole]

12880 DXDU=-AAA*SIN(U)*COS(V)*BZ(1,1)+BBB*COS(U)*COS(V)*BZ(1,2)

12890 DYDU=-AAA*SIN(U)*COS(V)*BZ(2,1)+BBB*COS(U)*COS(V)*BZ(2,2)

12900 DZDU=-AAA*SIN(U)*COS(V)*BZ(3,1)+BBB*COS(U)*COS(V)*BZ(3,2)

12910 DXDV=-AAA*COS(U)*SIN(V)*BZ(1,1)-BBB*SIN(U)*SIN(V)*BZ(1,2)+CCC*COS(V)*BZ(1,3)

12920 DYDV=-AAA*COS(U)*SIN(V)*BZ(2,1)-BBB*SIN(U)*SIN(V)*BZ(2,2)+CCC*COS(V)*BZ(2,3)

12930 DZDV=-AAA*COS(U)*SIN(V)*BZ(3,1)-BBB*SIN(U)*SIN(V)*BZ(3,2)+CCC*COS(V)*BZ(3,3)

12940 GOTO 13120

12950 '

12960 '[j-pole]

12970 DXDU=AAA*COS(U)*COS(V)*BZ(1,1)-CCC*SIN(U)*COS(V)*BZ(1,3)

12980 DYDU=AAA*COS(U)*COS(V)*BZ(2,1)-CCC*SIN(U)*COS(V)*BZ(2,3)

12990 DZDU=AAA*COS(U)*COS(V)*BZ(3,1)-CCC*SIN(U)*COS(V)*BZ(3,3)

13000 DXDV=-AAA*SIN(U)*SIN(V)*BZ(1,1)+BBB*COS(V)*BZ(1,2)-CCC*COS(U)*SIN(V)*BZ(1,3)

13010 DYDV=-AAA*SIN(U)*SIN(V)*BZ(2,1)+BBB*COS(V)*BZ(2,2)-CCC*COS(U)*SIN(V)*BZ(2,3)

13020 DZDV=-AAA*SIN(U)*SIN(V)*BZ(3,1)+BBB*COS(V)*BZ(3,2)-CCC*COS(U)*SIN(V)*BZ(3,3)

13030 GOTO 13120

13040 '

13050 '[i-pole]

13060 DXDU=-BBB*SIN(U)*COS(V)*BZ(1,2)+CCC*COS(U)*COS(V)*BZ(1,3)

13070 DYDU=-BBB*SIN(U)*COS(V)*BZ(2,2)+CCC*COS(U)*COS(V)*BZ(2,3)

13080 DZDU=-BBB*SIN(U)*COS(V)*BZ(3,2)+CCC*COS(U)*COS(V)*BZ(3,3)

13090 DXDV=AAA*COS(V)*BZ(1,1)-BBB*COS(U)*SIN(V)*BZ(1,2)-CCC*SIN(U)*SIN(V)*BZ(1,3)

13100 DYDV=AAA*COS(V)*BZ(2,1)-BBB*COS(U)*SIN(V)*BZ(2,2)-CCC*SIN(U)*SIN(V)*BZ(2,3)

13110 DZDV=AAA*COS(V)*BZ(3,1)-BBB*COS(U)*SIN(V)*BZ(3,2)-CCC*SIN(U)*SIN(V)*BZ(3,3)

13120 '

13130 DX=DYDU*DZDV-DYDV*DZDU

13140 DY=DZDU*DXDV-DZDV*DXDU

13150 DZ=DXDU*DYDV-DXDV*DYDU

13160 RETURN

13170 '

13180 ' *** 座標軸 ***

13190 '

13200 *DRAW.AXIS.WF:

13210 CLS 3

13220 SX1=59:SY1=360*(1-RY):SX2=580*RX+59:SY2=360

13230 LINE(SX1,SY1)-(SX2,SY2),GRP.COLOR,B

13240 RR=3:GOSUB *X.SCALE2:GOSUB *X.REF2

13250 RR=3:GOSUB *Y.SCALE2:GOSUB *Y.REF2

13260 I=IB:J=JB:RR=3:GOSUB *DOT.PLOT

13270 GOSUB *TITLE.BACK2.WF

13280 RETURN

13290 '

13300 ' *** タイトル ***

13310 '

13320 *TITLE.BACK2.WF:

13330 LX=INT((580*RX+59)/8)

13340 LY=INT(360*(1-RY)/16)

13350 LZ=INT((290*RX+59)/8)

13360 L1$="ellipsoid projection of "

13370 L2$="v"+MID$(STR$(IB),2)+"-"

13380 L3$="v"+MID$(STR$(JB),2)+"-"

13390 L4$="v"+MID$(STR$(KB),2)+" space"

13400 L$=L1$+L2$+L3$+L4$

13410 LOCATE LZ-LEN(L$)\2,LY-3:PRINT L$

13420 L$=UNIT$(JB)

13430 IF LEN(L$)=0 THEN 11090

13440 LOCATE 7,LY-1:PRINT "("+L$+")"

13450 L$=UNIT$(IB)

13460 IF LEN(L$)=0 THEN 11120

13470 LOCATE LX+3,21:PRINT "("+L$+")"

13480 '

13490 LOCATE LX+3,LY:PRINT "v components: "

13500 FOR K=1 TO M0

13510　 V$="v"+MID$(STR$(K),2)+"="

13520　 LOCATE LX+3,LY+K+1:PRINT V$;:PRINT V(K)

13530 NEXT K

13540 RETURN

13550 '

13560 ' *** 楕円体三面図 ***

13570 '

13580 *ELLIPSOID3:

13590 IF M<3 THEN RETURN

13600 '

13610 CLS 3

13620 IJK=2

13630 IB=VAL(MID$(KEY1$,1,1))

13640 JB=VAL(MID$(KEY1$,2,1))

13650 KB=VAL(MID$(KEY1$,3,1))

13660 BIAS=0

13670 IF IJK<>0 THEN GOSUB *EIGEN.WF

13680 IF IJK=0 AND STAGE$="2" THEN GOSUB *SECTION3

13690 IF IJK=0 THEN GOSUB *COMMON3:GOTO 13800

13700 '

13710 SX1=59:SY1=360*(1-RY/2):SX2=580*RX/2+59:SY2=360

13720 VHP=1:I=IB:J=JB:GOSUB *ELLIPSOID3.SUB

13730 '

13740 SX1=59:SY1=360*(1-RY):SX2=580*RX/2+59:SY2=360*(1-RY/2)

13750 VHP=2:I=IB:J=KB:GOSUB *ELLIPSOID3.SUB

13760 '

13770 SX1=580*RX/2+59:SY1=360*(1-RY/2):SX2=580*RX+59:SY2=360

13780 VHP=3:I=KB:J=JB:GOSUB *ELLIPSOID3.SUB

13790 'DRAW.COLOR=1:VHP=0:IJK=0:GOSUB *COMMON3:DRAW.COLOR=4

13800 '

13810 GOSUB *TITLE.BACK3:'FL$="b:e3":GOSUB *GSAVE

13820 GOSUB *TEMP

13830 RETURN

13840 '

13850 *COMMON3:

13860 SX1=59:SY1=360*(1-RY/2):SX2=580*RX/2+59:SY2=360

13870 I=IB:J=JB:GOSUB *ELLIPSOID3.SUB

13880 '

13890 SX1=59:SY1=360*(1-RY):SX2=580*RX/2+59:SY2=360*(1-RY/2)

13900 I=IB:J=KB:GOSUB *ELLIPSOID3.SUB

13910 '

13920 SX1=580*RX/2+59:SY1=360*(1-RY/2):SX2=580*RX+59:SY2=360

13930 I=KB:J=JB:GOSUB *ELLIPSOID3.SUB

13940 'SX1=580*RX/2+59:SY1=360*(1-RY):SX2=580*RX+59:SY2=360*(1-RY/2)

13950 RETURN

13960 '

13970 *ELLIPSOID3.SUB:

13980 UO=VO(IB):UL=VL(IB)

13990 VO=VO(JB):VL=VL(JB)

14000 DU=VD(IB):SU=VS(IB)

14010 DV=VD(JB):SV=VS(JB)

14020 '

14030 LINE(SX1,SY1)-(SX2,SY2),GRP.COLOR,B

14040 RR=2:GOSUB *X.SCALE2

14050 RR=2:GOSUB *Y.SCALE2

14060 'RR=2:GOSUB *DOT.PLOT

14070 GOSUB *DRAW.ELLIPSE3

14080 RETURN

14090 '

14100 ' *** 三面鏡 ***

14110 '

14120 *DRAW.ELLIPSE3:

14130 'DRAW.COLOR=4

14140 WUO=UO-UL:WUL=UO+UL

14150 WVO=VO-VL:WVL=VO+VL

14160 WINDOW(WUO,-WVL)-(WUL,-WVO)

14170 VIEW(SX1,SY1)-(SX2,SY2)

14180 '

14190 IF IJK<>0 THEN GOSUB *COMMON.WF

14200 IF IJK=0 AND STAGE$="1" THEN GOSUB *CONTOUR

14210 IF IJK=0 AND STAGE$="2" THEN GOSUB *CONTOUR

14220 'IF IJK=0 AND STAGE$="2" THEN GOSUB *SECTION

14230 'GOSUB *COMMON.WF

14240 'GOSUB *CONTOUR

14250 'GOSUB *SECTION

14260 '

14270 WINDOW(0,0)-(639,399)

14280 VIEW(0,0)-(639,399)

14290 RETURN

14300 '

14310 ' *** ３次元切断面の輪郭 ***

14320 '

14330 *SECTION3:

14340 AZ(1,1)=HESSE(IB,IB)

14350 AZ(1,2)=HESSE(IB,JB)

14360 AZ(1,3)=HESSE(IB,KB)

14370 AZ(2,1)=AZ(1,2)

14380 AZ(2,2)=HESSE(JB,JB)

14390 AZ(2,3)=HESSE(JB,KB)

14400 AZ(3,1)=AZ(1,3)

14410 AZ(3,2)=AZ(2,3)

14420 AZ(3,3)=HESSE(KB,KB)

14430 NZ=3:GOSUB *INVERSE

14440 '

14450 BIAS=0:' [*]

14460 FOR I=1 TO 3

14470　 FOR J=1 TO 3:COV(I+BIAS,J+BIAS)=BZ(I,J):NEXT J

14480 NEXT I

14490 RETURN

14500 '

14510 ' *** タイトル ***

14520 '

14530 *TITLE.BACK3:

14540 LX=INT((580*RX+59)/8)

14550 LY=INT(360*(1-RY)/16)

14560 LZ=INT((290*RX+59)/8)

14570 L1$="ellipsoid projection of "

14580 L2$="v"+MID$(STR$(IB),2)+"-"

14590 L3$="v"+MID$(STR$(JB),2)+"-"

14600 L4$="v"+MID$(STR$(KB),2)+" space"

14610 L$=L1$+L2$+L3$+L4$

14620 LOCATE LZ-LEN(L$)\2,LY-3:PRINT L$

14630 '

14640 LOCATE LX+3,LY:PRINT "v components: "

14650 FOR K=1 TO M0

14660　 V$="v"+MID$(STR$(K),2)+"="

14670　 LOCATE LX+3,LY+K+1:PRINT V$;:PRINT V(K)

14680 NEXT K

14690 '

14700 LOCATE LX+3,19:PRINT "horizontal|"

14710 LOCATE LX+3,20:PRINT "-------------------"

14720 LOCATE LX+3,21:PRINT " vertical | profile"

14730 RETURN

14740 '

14750 ' ** VRAM SAVE **

14760 '

14770 *GSAVE:

14780 DEF SEG=&HA000

14790 BSAVE FL$+".CHR",&H0,&HFA0

14800 DEF SEG=&HA200

14810 BSAVE FL$+".ATR",&H0,&HFA0

14820 DEF SEG=&HA800

14830 BSAVE FL$+".BLU",&H0,&H7D00

14840 DEF SEG=&HB000

14850 BSAVE FL$+".RED",&H0,&H7D00

14860 DEF SEG=&HB800

14870 BSAVE FL$+".GRN",&H0,&H7D00

14880 RETURN

　

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【３】プログラムの概説

　

　コラム「高次元の正多胞体」が「ロボットアームと６次元直方体」の続編だったように，このシリーズは「ロボットアームと６次元楕円体」の続編です．

　

　前編「ロボットアームと６次元楕円体」で紹介したプログラムは，

　　(1)６次元楕円を直接指定した２次元平面（３次元空間）へ投影する

ものだったのですが，続編ではそれを

　　(a)６次元楕円の任意の超平面での切り口をｍ次元空間に投影する（あるいは，ｍ次元空間での切り口か正射影を求める）

　　(b)それを指定した２次元平面（３次元空間）へ投影する

といった２ステップ方式の投影方法に変更しています．

　

　(a)を行う際，行列Ｍを，ｍ次元空間での切り口か正射影かに応じて

　　α2~α1＝０ 　　→　Ｍ＝［Ａ］

　　α2~α1＝Ｄ~Ｃ　→　Ｍ＝［Ａ−ＢＤ~Ｃ］

とすれば，切り口の正射影を求めたり，正射影の切り口を求めたりすることもできるようになります．

　

　６次元楕円の２次元平面および３次元空間への正射影を考えると，３次元空間への正射影（３次元楕円）を２次元平面に投影した楕円と，最初から２次元平面に正射影した楕円は一致します．

　

　それに対して，６次元楕円の２次元平面での切り口と３次元空間での切り口（３次元楕円）を２次元平面に投影したものでは，楕円の大きさが一致しません．２次元切り口は３次元切り口の２次元平面での切り口であって，３次元切り口の２次元平面への正射影ではないからです．

　

　歯切れの悪い説明で申し訳ありませんが，正射影の正射影は正射影ですが，切り口の切り口は切り口であって正射影ではないために，楕円の大きさが一致しないのです．

　

　プログラム（*hyperplane）では，ｍ次元空間の定義を

　　gosub *hyperplane.def

で行っています．それに従って，

　　gosub *hyperplane.1

で，分割型正方行列

　　Ｈ＝［Ａ，Ｂ］

　　　　［Ｃ，Ｄ］

を取り込んで，

　　gosub *hyperplane.2,

　　gosub *hyperplane.3

で，任意の超平面による切り口のｍ次元空間への投影：

　　Ｍ＝［Ａ−Ｂα2~α1−（α2~α1）’Ｃ＋（α2~α1）’Ｄα2~α1］

を計算しています．

　

　また，gosub *hyperplaneの代わりに

　　gosub *hyper.contour

とすればｍ次元空間への正射影：Ｍ＝［Ａ−ＢＤ~Ｃ］，

　　gosub *hyper.section

ではｍ次元空間での切り口：Ｍ＝［Ａ］となります．

　

　これらの考え方については（その４）で既に説明してありますが，実際のプログラム作りに際しては，(b)についての書き直しが最小になるように配慮しました．

　

　実際に動かしてみると，*hyperplaneは思ったほどには機能せず，アイディア倒れの感なきにしもあらずでしたが，*hyper.contourや*hyper.sectionは以前作成したプログラムとの比較を行う上で，有益な情報をもたらしてくれました．

　

　また，「ロボットアームと６次元楕円体」に掲げたプログラムでは，６次元楕円を直接２次元平面に投影するため，３次元楕円の極と元の６次元楕円の極は無関係でしたが，いったん３次元空間に投影することによって，２次元楕円上に３次元楕円の極を描くことができるようになりました．

　

　なお，このプログラムは９次元以下であれば修正なしに使えますが，１０次元を超えると文字列処理の部分を書き直す必要が生じます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝