シリコンフラーレンの木工模型は切頂ねじれ重角錐に近似することによってできると思われましたが，川添良幸先生から頂いた数値だけでは曲がっている五辺形面を平面でうまく近似することはできず，そのため，（その１）では上下の四角錐台（蓋と底）と中央の四角反柱（胴）の三段に分けて作るしかありませんでした．
３ブロック模型では不満が残る結果となってしまったわけですが，ｎ^2５^2nは重要ですから，どうしても４^2５^8の１ブロック模型を作りたかったのです．これを解決するには曲面五辺形を近似する平面五角形を求めなければなりません．
しかし，単純な最小２乗法による平面近似では閉じた形のねじれ重角錐は組めませんので，閉じた形に組めるための条件を満足させつつ最小２乗法を適用することになります．また，１面だけの近似では隣接する面にも近似の影響が及んできますから，１０面すべてを同時に４^2５^8型にフィットさせる必要もでてきます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】空間の回転
最小２乗法による平面近似のためには頂点の（ｘ，ｙ，ｚ）座標が必要になります．クマール教授に送っていただいたシリコンフラーレンの座標データを掲げます．
Dear Dr. Sato,
Thank you very much. I am happy that you are helping to prepare a wood model for us. Yes I mentioned this to Prof Kawazoe that the pentagons are not planar.
Here are the coordinates and I hope this would help you. Let me know if you could work with these.
Sincerely yours
Vijay Kumar
1.934516430 1.542401552 1.346459150
1.019369125 -0.1833969355 2.539282084
1.066837907 1.645336747 -2.071133852
-2.042828560 2.019238472 1.917411089
-0.9660972357 -1.590861082 -1.912511110
-1.671085477 -1.770418882 1.542118192
1.214756370 -1.702803135 -0.9662210941
-2.162976742 2.027619362 -1.608949900
-1.260006189 0.8881688118E-01 2.866134405
-0.8372278214 0.5248215795 -2.775815010
0.5208492279E-01 2.999498367 1.355939865
0.6083330512 -2.041862488 1.214753628
-0.2587424219 3.149661779 -0.9066201448
2.356674671 0.3653883934 -0.6785058975
-3.271762133 1.036553144 0.1317510009
-2.763231754 -1.227902293 -0.3948068619
-0.4362684488 0.4303944111 0.1003119349
これらを（ｘi，ｙi，ｚi）としますが，このうちi=17はシリコンフラーレン内部に閉じこめられた金属元素と思われます．i=1~16がシリコンの座標データですが，このままでは扱いにくいので，重心を中心として上蓋が正方形，下蓋がそれを４５°回転させた標準的な座標系に変換する必要があります．
そのためには，まず点Ｐ（ｘ，ｙ，ｚ）をｚ軸の周りにθ，ｘ軸の周りにφ回転させることを考えます．
（１，０，０）→（ｃｏｓθ，ｓｉｎθ，０）
→（ｃｏｓθ，ｓｉｎθｃｏｓφ，ｓｉｎθｓｉｎφ）
（０，１，０）→（−ｓｉｎθ，ｃｏｓθ，０）
→（−ｓｉｎθ，ｃｏｓθｃｏｓφ，ｃｏｓθｓｉｎφ）
（０，０，１）→（０，０，１）
→（０，−ｓｉｎφ，ｃｏｓφ）
ですから，Ｐ（ｘ，ｙ，ｚ）は
［ｃｏｓθ ，−ｓｉｎθ ， ０］［ｘ］
［ｓｉｎθｃｏｓφ，ｃｏｓθｃｏｓφ，−ｓｉｎφ］［ｙ］
［ｓｉｎθｓｉｎφ，ｃｏｓθｓｉｎφ， ｃｏｓφ］［ｚ］
に移ります．
行列
［ｃｏｓθ ，−ｓｉｎθ ， ０］
［ｓｉｎθｃｏｓφ，ｃｏｓθｃｏｓφ，−ｓｉｎφ］
［ｓｉｎθｓｉｎφ，ｃｏｓθｓｉｎφ， ｃｏｓφ］
は３次元の直交行列です．→［補］
次に，ベクトルＯＰ↑＝（ｘ，ｙ，ｚ）を新しいｘ軸とします．ｘ^2＋ｙ^2＋ｚ^2＝ａ^2とすると方向余弦は（ｘ／ａ，ｙ／ａ，ｚ／ａ）で，このとき
ｃｏｓθ＝ｘ／ａ，ｓｉｎθ＝√（ｙ^2＋ｚ^2）／ａ
また，ｔａｎφ＝ｚ／ｙですから
ｃｏｓφ＝ｙ／√（ｙ^2＋ｚ^2），ｓｉｎφ＝ｙ／√（ｙ^2＋ｚ^2）
なる幾何学的関係が成り立ちます（象限判別が必要となるが割愛）．
空間の回転により，旧座標系の点（ａ，０，０）が点（ｘ，ｙ，ｚ）すなわち点（１，０，０）が点（ｘ／ａ，ｙ／ａ，ｚ／ａ）に移ることになりますから，
［ｘ］ ［ｃｏｓθ ，−ｓｉｎθ ， ０］［ａ］
［ｙ］＝［ｓｉｎθｃｏｓφ，ｃｏｓθｃｏｓφ，−ｓｉｎφ］［０］
［ｚ］ ［ｓｉｎθｓｉｎφ，ｃｏｓθｓｉｎφ， ｃｏｓφ］［０］
です．
この点が新しい座標系における点［ａ，０，０］になるわけですから，この直交行列の逆行列をとると，旧座標系の（ｘ，ｙ，ｚ）と新座標系の［Ｘ，Ｙ，Ｚ］の間には
［Ｘ］ ［ ｃｏｓθ，ｓｉｎθｃｏｓφ，ｓｉｎθｓｉｎφ］［ｘ］
［Ｙ］＝［−ｓｉｎθ，ｃｏｓθｃｏｓφ，ｃｏｓθｓｉｎφ］［ｙ］
［Ｚ］ ［ ０， ｓｉｎφ， ｃｏｓφ］［ｚ］
なる変換式が成立することになります．
こうして得られた点をあらためてＰ（ｘ，ｙ，ｚ）とします．実際に変換された点の座標をみますと，天地面は真の正方形ではないことがわかります．そうなると五角形面も１種類ではないことになりますが，些細な点には目をつぶることにします．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】４^2５^8型の目的関数
（その２）で述べた噛み合わせ可能な五角形の頂点をＡとすると，辺ＣＤが正ｎ角形と組み合わさる辺でその長さをａ，辺ＢＣ＝辺ＤＥ＝ｂ，辺ＡＢ＝辺ＥＡ＝ｃ，底角∠Ｃ＝∠Ｄ＝θ，頂角∠Ａ＝φで表します．θは９０°＜θ＜（９０＋１８０／ｎ）°の値をとるものとします．
正ｎ角形面の中心から辺までの距離ｄは
ｄ＝ａ／（２ｓｉｎ(π／ｎ)）
また，
ｙ＝−ｃｏｓθ／ｓｉｎ(π／ｎ)
ｚ＝−（ｂ^2−ｙ^2）^(1/2)
ζ＝z/y{(y+a/(2sin(π/n))(tan(π/n)sin(π/n)+cos(π/n)-1)}+z
とおくと，正ｎ角形面間の距離ｈは
ｈ＝｜ｚ＋ζ｜
で与えられます．
たとえば，頂点Ｃを（−ａ／，ａ／２，ｈ／２），頂点Ｄを（ａ／２，ａ／２，ｈ／２）にとると，
Ａ（０，ｙ＋ｄ，ｈ／２＋ζ）
Ｂ（−（ｙ＋ｄ）ｓｉｎ(π／ｎ)，（ｙ＋ｄ）ｃｏｓ(π／ｎ)，ｈ／２＋ｚ）
Ｄ（（ｙ＋ｄ）ｓｉｎ(π／ｎ)，（ｙ＋ｄ）ｃｏｓ(π／ｎ)，ｈ／２＋ｚ）
となります．以下同様の要領でシリコンフラーレンの１６点すべてをａ，ｂ，θでパラメトライズして得られる点をＱ（ｘ，ｙ，ｚ）とします．
多面体の空間データの場合はｘ軸，ｙ軸，ｚ軸のいずれにも誤差がありますから，目的関数を
ｓ＝Σ（Ｑi−Ｐi）^2＝Σ｛(qx-px)^2+(qy-py)^2+(qz-pz)^2｝
として最小２乗法を適用すると１０面すべてを同時に４^2５^8型フィットさせることができます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】計算例
（その２）においてシリコンフラーレン型（４^2５^8）では等辺（ａ＝ｂ＝ｃ）のものはありえないことが証明されています．そこで，可変パラメータとしてａ，ｂ，θを選び，五角形面が噛み合うための条件からｃ，φを決定しつつ最小２乗法を適用することになります．
実際のプログラムでは微分行列（ヤコビアン）を用いないで直接探索するシンプレックス法で最小２乗近似しました．また，本来は負になりえないパラメータが負の値となることを避けるために，パラメータに上限・下限を設けてパラメータが常にその範囲内に収まるような制約条件付き最小２乗近似を行っています．
ａ＝１に換算した計算結果を示しますと
ａ＝１，ｂ＝．９６０８５１，ｃ＝１．２０６２９，ｈ＝２．１１６１９
θ＝１１７．１９，φ＝１０２．２４
となりました．このとき二面角は１２０．９１１と計算されます．また，ａ＝ｂなる条件下では
ａ＝１，ｂ＝１，ｃ＝１．２２９０８，ｈ＝２．１７１２１
θ＝１１７．４５２，φ＝１０２．８６８，二面角＝１２１．２９９
となりました．
いずれの場合も二面角は１２１°前後，また，底角は１１７．３°前後ということになり，微妙な差しかありませんでした．それに対して，川添先生から頂いたデータでは底角は１１２．２°でしたから約５°の違いがあります．以下に中川宏さん作によるシリコンフラーレンの１ブロック模型を掲げます．
３ブロック模型と較べてまあまあよく似ているという印象です．最小２乗法を用いて最適なものを選んだとはいうものの，極端な曲面であればそもそも平面近似自体が無理なわけです．幸い，シリコンフラーレンでは角不足はわずか３°であり，そのためにうまく近似できたのだと思われます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】プログラム
1000 '
1010 ' **** silicon fullerene (SF) ****
1020 ' 2005/11/05 (C) ｻﾄｳ ｲｸﾛｳ
1030 ON STOP GOSUB *STOP.RTN:STOP ON
1040 SCREEN 3,0:CONSOLE ,,0,1
1050 CLS 3:WIDTH 80,25:COLOR 6
1060 GOSUB *INITIALIZE
1070 GOSUB *READ.DATA
1080 GOSUB *SIMPLEX
1090 END
1100 '
1110 *STOP.RTN:
1120 CLOSE:SCREEN 3,0:CONSOLE ,,0,1:CLS 3
1130 STOP
1140 '
1150 ' *** 初期設定 ***
1160 '
1170 *INITIALIZE:
1180 PI=3.14159
1190 '
1200 N=4
1210 M=3
1220 N1=N*4
1230 '
1240 DIM X(N1),Y(N1),Z(N1)
1250 DIM XX(N1),YY(N1),ZZ(N1)
1260 DIM P(M),PL(M),PU(M)
1270 DIM Q(M),CN(M)
1280 DIM VX(M+1,M),SS(M+1)
1290 DIM VR(M),VC(M),VE(M)
1300 RETURN
1310 '
1320 ' *** データの読み込み ***
1330 '
1340 *READ.DATA:
1350 RESTORE *J
1360 FOR I=1 TO N1
1370 READ X(I),Y(I),Z(I)
1380 NEXT I
1390 X=0:Y=0:Z=0
1400 FOR I=1 TO N1
1410 X=X+X(I)
1420 Y=Y+Y(I)
1430 Z=Z+Z(I)
1440 NEXT I
1450 X0=X/N1:Y0=Y/N1:Z0=Z/N1
1460 '
1470 FOR I=1 TO N1
1480 X(I)=X(I)-X0
1490 Y(I)=Y(I)-Y0
1500 Z(I)=Z(I)-Z0
1510 NEXT I
1520 '
1530 X=X(2)-X(1)
1540 Y=Y(2)-Y(1)
1550 Z=Z(2)-Z(1)
1560 AA=SQR(X*X+Y*Y+Z*Z)
1570 YZ=SQR(Y*Y+Z*Z)
1580 C1=X/AA
1590 S1=YZ/AA
1600 C2=Y/YZ
1610 S2=Z/YZ
1620 '
1630 FOR I=1 TO N1
1640 X=X(I):Y=Y(I):Z=Z(I)
1650 X(I)=C1*X+S1*C2*Y+S1*S2*Z
1660 Y(I)=-S1*X+C1*C2*Y+C1*S2*Z
1670 Z(I)=-S2*Y+C2*Z
1680 NEXT I
1690 FOR I=1 TO N1
1700 PRINT X(I),Y(I),Z(I)
1710 NEXT I
1720 '
1730 P(1)=AA
1740 P(2)=AA
1750 P(3)=(1+1/N)*PI/2
1760 '
1770 CN(1)=3:CN(2)=3:CN(3)=3
1780 PL(1)=P(1)*.8:PU(1)=P(1)*1.2
1790 PL(2)=P(2)*.8:PU(2)=P(2)*1.2
1800 PL(3)=PI/2 :PU(3)=(1/2+1/N)*PI
1810 RETURN
1820 '
1830 *J:
1840 DATA -.2587424219,3.149661779,-.9066201448 :'[13]
1850 DATA -2.162976742,2.027619362,-1.608949900 :'[8]
1860 DATA -.8372278214,.5248215795,-2.775815010 :'[10]
1870 DATA 1.066837907,1.645336747,-2.071133852 :'[3]
1880 '
1890 DATA -2.042828560,2.019238472,1.917411089 :'[4]
1900 DATA -2.763231754,-1.227902293,-.3948068619 :'[16]
1910 DATA 1.214756370,-1.702803135,-.9662210941 :'[7]
1920 DATA 1.934516430,1.542401552,1.346459150 :'[1]
1930 '
1940 DATA -3.271762133,1.036553144,.1317510009 :'[15]
1950 DATA -.9660972357,-1.590861082,-1.912511110 :'[5]
1960 DATA 2.356674671,.3653883934,-.6785058975 :'[14]
1970 DATA .5208492279E-01,2.999498367,1.355939865 :'[11]
1980 '
1990 DATA -1.260006189,.8881688118E-01,2.866134405:'[9]
2000 DATA -1.671085477,-1.770418882,1.542118192 :'[6]
2010 DATA .6083330512,-2.041862488,1.214753628 :'[12]
2020 DATA 1.019369125,-.1833969355,2.539282084 :'[2]
2030 '
2040 ' ** PENTAGON **
2050 '
2060 *PENTAGON:
2070 AA=P(1)
2080 BB=P(2)
2090 C1=COS(P(3)):S1=SIN(P(3))
2100 '
2110 D=AA/2/SIN(PI/N)
2120 Y=-BB*C1/SIN(PI/N)
2130 Z=-SQR(BB*BB-Y*Y)
2140 C=(Y+D)*(Y+D)*((SIN(PI/N))^2+(1/COS(PI/N)-1)^2+(TAN(PI/N)*SIN(PI/N)+COS(PI/N)-1)^2*Z*Z/Y/Y)
2150 CC=SQR(C)
2160 C=1-2*(Y+D)*(Y+D)*SIN(PI/N)*SIN(PI/N)/CC/CC
2170 TH=ATN(SQR(1-C*C)/C)
2180 IF C<0 THEN TH=PI+TH
2190 '
2200 ZZ=Z/Y*(Y+D)*(TAN(PI/N)*SIN(PI/N)+COS(PI/N)-1)+Z
2210 HH=ABS(Z+ZZ)
2220 '
2230 XX(1)=-AA/2
2240 YY(1)=AA/2
2250 ZZ(1)=HH/2
2260 FOR K=1 TO N-1
2270 XX=XX(K):YY=YY(K)
2280 XX(1+K)= XX*COS(2*PI/N)+YY*SIN(2*PI/N)
2290 YY(1+K)=-XX*SIN(2*PI/N)+YY*COS(2*PI/N)
2300 ZZ(1+K)=HH/2
2310 NEXT K
2320 '
2330 XX(N+1)=0
2340 YY(N+1)=Y+D
2350 ZZ(N+1)=HH/2+ZZ
2360 FOR K=1 TO N-1
2370 XX=XX(N+K):YY=YY(N+K)
2380 XX(N+1+K)= XX*COS(2*PI/N)+YY*SIN(2*PI/N)
2390 YY(N+1+K)=-XX*SIN(2*PI/N)+YY*COS(2*PI/N)
2400 ZZ(N+1+K)=HH/2+ZZ
2410 NEXT K
2420 '
2430 XX(N*2+1)=YY(N+1)*SIN(PI/N)
2440 YY(N*2+1)=YY(N+1)*COS(PI/N)
2450 ZZ(N*2+1)=HH/2+Z
2460 FOR K=1 TO N-1
2470 XX=XX(N*2+K):YY=YY(N*2+K)
2480 XX(N*2+1+K)= XX*COS(2*PI/N)+YY*SIN(2*PI/N)
2490 YY(N*2+1+K)=-XX*SIN(2*PI/N)+YY*COS(2*PI/N)
2500 ZZ(N*2+1+K)=HH/2+Z
2510 NEXT K
2520 '
2530 XX(N*3+1)=0
2540 YY(N*3+1)=D
2550 ZZ(N*3+1)=-HH/2
2560 FOR K=1 TO N-1
2570 XX=XX(N*3+K):YY=YY(N*3+K)
2580 XX(N*3+1+K)= XX*COS(2*PI/N)+YY*SIN(2*PI/N)
2590 YY(N*3+1+K)=-XX*SIN(2*PI/N)+YY*COS(2*PI/N)
2600 ZZ(N*3+1+K)=-HH/2
2610 NEXT K
2620 RETURN
2630 '
2640 ' ** SUM of RESIDUALS **
2650 '
2660 *RESIDUAL:
2670 GOSUB *PENTAGON
2680 SS=0
2690 FOR K=1 TO N1
2700 SS=SS+(XX(K)-X(K))^2+(YY(K)-Y(K))^2+(ZZ(K)-Z(K))^2
2710 NEXT K
2720 RETURN
2730 '
2740 ' *** FINAL OUTPUT ***
2750 '
2760 *FINAL:
2770 SN=SS(IL)
2780 CLS 3
2790 PRINT "LOOP=";LOOP
2800 PRINT
2810 PRINT "PARAMETERS: "
2820 FOR I=1 TO M
2830 PRINT "P(";I;")=";P(I)
2840 NEXT I
2850 PRINT
2860 PRINT "RESIDUALS=";SN
2870 PRINT
2880 FOR I=1 TO N1
2890 PRINT XX(I),YY(I),ZZ(I)
2900 NEXT I
2910 '
2920 PRINT
2930 PRINT "何かキーを押してください"
2940 WHILE INKEY$="":WEND
2950 RETURN
2960 '
2970 ' *** シンプレックス法 ***
2980 '
2990 *SIMPLEX:
3000 LMAX=1000
3010 LOOP=0
3020 '
3030 VAA=1 :'[reflection]
3040 VBB=.5:'[contraction]
3050 VCC=2 :'[expansion]
3060 GOSUB *PARAM.SET
3070 *SIMPLEX.LOOP:
3080 GOSUB *SIMPLEX.EXECUTION
3090 GOSUB *SIMPLEX.TEMPORARY
3100 GOSUB *SIMPLEX.CONVERGENCE
3110 IF SW=1 AND LOOP<LMAX THEN GOTO *SIMPLEX.LOOP
3120 GOSUB *FINAL
3130 RETURN
3140 '
3150 ' *** PARAMETER SET ***
3160 '
3170 *PARAM.SET:
3180 GOSUB *VAR.CONV
3190 FOR J=1 TO M:VX(1,J)=Q(J):NEXT J
3200 FOR I=2 TO M+1
3210 FOR J=1 TO M
3220 VX(I,J)=(RND(1)*.1+1)*VX(1,J)
3230 NEXT J
3240 NEXT I
3250 '
3260 FOR I=1 TO M+1
3270 FOR J=1 TO M
3280 Q(J)=VX(I,J)
3290 NEXT J
3300 GOSUB *REV.CONV
3310 GOSUB *RESIDUAL
3320 SS(I)=SS
3330 NEXT I
3340 RETURN
3350 '
3360 ' *** EXECUTION of FITTING ***
3370 '
3380 *SIMPLEX.EXECUTION:
3390 IH=1:IL=1
3400 FOR I=2 TO M+1
3410 IF SS(I)>SS(IH) THEN IH=I
3420 IF SS(I)<SS(IL) THEN IL=I
3430 NEXT I
3440 SH=SS(IH):SL=SS(IL)
3450 FOR J=1 TO M
3460 VSUM=0
3470 FOR I=1 TO M+1
3480 IF I=IH THEN 3500
3490 VSUM=VSUM+VX(I,J)
3500 NEXT I
3510 VG=VSUM/M
3520 VR(J)=(1+VAA)*VG-VAA*VX(IH,J)
3530 VC(J)=VBB*VX(IH,J)+(1-VBB)*VG
3540 VE(J)=VCC*VR(J)+(1-VCC)*VG
3550 NEXT J
3560 FOR J=1 TO M:Q(J)=VR(J):NEXT J
3570 GOSUB *REV.CONV
3580 GOSUB *RESIDUAL:SR=SS
3590 IF SR<SL THEN GOSUB *CASE1 ELSE IF SR<SH THEN GOSUB *CASE2 ELSE GOSUB *CASE3
3600 RETURN
3610 '
3620 *CASE1:
3630 FOR J=1 TO M:Q(J)=VE(J):NEXT J
3640 GOSUB *REV.CONV
3650 GOSUB *RESIDUAL:SE=SS
3660 IF SE<SR THEN GOSUB *REPLACE:GOTO 3690
3670 FOR J=1 TO M:VX(IH,J)=VR(J):NEXT J
3680 SS(IH)=SR
3690 RETURN
3700 '
3710 *CASE2:
3720 GOSUB *REPLACE
3730 RETURN
3740 '
3750 *CASE3:
3760 FOR J=1 TO M:Q(J)=VC(J):NEXT J
3770 GOSUB *REV.CONV
3780 GOSUB *RESIDUAL:SC=SS
3790 IF SC<SH THEN GOSUB *REPLACE:GOTO 3910
3800 FOR I=1 TO M+1
3810 FOR J=1 TO M:Q(J)=VX(IL,J):NEXT J
3820 IF I=IL THEN 3900
3830 FOR J=1 TO M
3840 VX(I,J)=.5*(VX(I,J)+VX(IL,J))
3850 Q(J)=VX(I,J)
3860 NEXT J
3870 GOSUB *REV.CONV
3880 GOSUB *RESIDUAL
3890 SS(I)=SS
3900 NEXT I
3910 RETURN
3920 '
3930 *REPLACE:
3940 FOR J=1 TO M
3950 VX(IH,J)=Q(J)
3960 NEXT J
3970 SS(IH)=SS
3980 RETURN
3990 '
4000 ' *** TEMPORARY OUTPUT ***
4010 '
4020 *SIMPLEX.TEMPORARY:
4030 FOR J=1 TO M:Q(J)=VX(IL,J):NEXT J
4040 SS=SS(IL)
4050 GOSUB *REV.CONV
4060 '
4070 LOOP=LOOP+1
4080 CLS 3
4090 PRINT "LOOP=";LOOP
4100 PRINT
4110 PRINT "PARAMETERS: "
4120 FOR J=1 TO M
4130 PRINT "P(";J;")=";P(J)
4140 NEXT J
4150 PRINT
4160 PRINT "RESIDUALS=";SS
4170 RETURN
4180 '
4190 ' *** CHECK of CONVERGENCE ***
4200 '
4210 *SIMPLEX.CONVERGENCE:
4220 EPS=.0001
4230 SW=0
4240 FOR J=1 TO M
4250 VMAX=VX(1,J):VMIN=VMAX
4260 FOR I=2 TO M+1
4270 IF VX(I,J)>VMAX THEN VMAX=VX(I,J)
4280 IF VX(I,J)<VMIN THEN VMIN=VX(I,J)
4290 NEXT I
4300 IF ABS((VMAX-VMIN)/VMIN)>EPS THEN SW=1
4310 NEXT J
4320 RETURN
4330 '
4340 ' *** 変数変換 ***
4350 '
4360 *VAR.CONV:
4370 FOR IS=1 TO M
4380 CN=CN(IS):PM=P(IS):PL=PL(IS):PU=PU(IS)
4390 IF CN=0 THEN Q(IS)=PM
4400 IF CN=1 THEN Q(IS)=SQR(PM)
4410 IF CN=2 THEN Q(IS)=LOG(PM)
4420 IF CN=3 THEN TM=SQR((PM-PL)/(PU-PL)):Q(IS)=ATN(TM/SQR(1-TM*TM))
4430 IF CN=4 THEN Q(IS)=LOG((PM-PL)/(PU-PM))/2
4440 NEXT IS
4450 RETURN
4460 '
4470 ' *** 変数逆変換 ***
4480 '
4490 *REV.CONV:
4500 FOR IS=1 TO M
4510 CN=CN(IS):QM=Q(IS):PL=PL(IS):PU=PU(IS)
4520 IF CN=0 THEN P(IS)=QM
4530 IF CN=1 THEN P(IS)=QM*QM
4540 IF CN=2 THEN P(IS)=EXP(QM)
4550 IF CN=3 THEN P(IS)=PL+(PU-PL)*SIN(QM)*SIN(QM)
4560 IF CN=4 THEN P(IS)=PL+(PU-PL)*EXP(QM)/(EXP(QM)+EXP(-QM))
4570 NEXT IS
4580 RETURN
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【補】直交行列と回転行列
二次元の場合，直交行列は必ず原点の周りの回転
R=[cosθ,-sinθ]
[sinθ, cosθ]
で表現可能であるので，平面の回転では回転行列Ｒをかけてやればよいのですが，空間の回転では直交行列は回転行列とは限りません．
たとえば，
[1 0 0]
[0 1 0]
[0 0 -1]
はユニタリではあるが，軸の周りの「回転」ではありません．（教訓）：２次元の直交変換はそのまま平面の回転なのですが，３次元の直交変換は一般に空間の回転になりません．
そこで，空間を回転させる行列で直交変換となっているものを探してみたところ，パラメータ数が３つの「回転」かつ「直交」行列として
(1)オイラー角に基づくもの
(2)ロール・ピッチ・ヨーに基づくもの
が見つかりました．
(1)はｚ軸まわりの回転α→新しいｙ軸まわりの回転β→新しいｚ軸まわりの回転γ，(2)はｚ軸まわりの回転φ→新しいｙ軸まわりの回転θ→新しいｘ軸まわりの回転ψの３段階によって表すもので，３番目の回転軸が異なるだけで両者に本質的な違いはありません．
ロボット工学の基礎としては(2)が良く扱われるとのことですから，ここでは(2)を示しますが，
Ｒψ＝［１，０，０］
［０，ｃｏｓψ，−ｓｉｎψ］
［０，ｓｉｎψ，ｃｏｓψ］
Ｒθ＝［ｃｏｓθ，０，−ｓｉｎθ］
［０，１，０］
［ｓｉｎθ，０，ｃｏｓθ］
Ｒφ＝［ｃｏｓφ，−ｓｉｎφ，０］
［ｓｉｎφ，ｃｏｓφ，０］
［０，０，１］
の合成マトリックス：Ｒ＝ＲφＲθＲψは
R(1,1)=cosφcosθ
R(2,1)=sinφcosθ
R(3,1)=-sinθ
R(1,2)=cosφsinθsinψ-sinφcosψ
R(2,2)=sinφsinθsinψ+cosφcosψ
R(3,2)=cosθsinψ
R(1,3)=cosφsinθcosψ+sinφsinψ
R(2,3)=sinφsinθcosψ-cosφsinψ
R(3,3)=cosθcosψ
のようになります．
R=[cosθ,-sinθ]
[sinθ, cosθ]
を一般化したものですが，非常に面倒な表現になってしまいました．それでも回転行列を与えられた場合のφ，θ，ψの解法は知られているようです（ただし，一意ではなく，象限判別が必要となる）．
また，単位ベクトル
ｎ↑＝（α,β,γ） α，β，γは方向余弦で，α^2+β^2+γ^2＝１を満たすものとする
を回転軸とし，その周りに正の回転方向にθだけ回転する回転行列
R(1,1)=α^2(1-cosθ)+cosθ
R(2,2)=β^2(1-cosθ)+cosθ
R(3,3)=γ^2(1-cosθ)+cosθ
R(1,2)=αβ(1-cosθ)+γsinθ
R(2,1)=αβ(1-cosθ)-γsinθ
R(1,3)=αγ(1-cosθ)-βsinθ
R(3,1)=αγ(1-cosθ)+βsinθ
R(2,3)=βγ(1-cosθ)+αsinθ
R(3,2)=βγ(1-cosθ)-αsinθ
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
