３．3Dの直線と平面

3Dの直線と平面

1. 直線

2Dであれば傾きと切片が直線を定め、xに対するyの値が計算できます。

3Dであっても直線が2点で決まることは納得できます。しかし、何に対して何を導くのか良くわかりません。

どうやら、直線状の点 $P({x}_{p},{y}_{p},{z}_{p})$ を媒介変数tで表すようです。2Dの直線は1つの式で表されますが、3Dの直線は3つの式で表すようです。

点Aを通る方向ベクトルdの直線上の点Pは、A、Pをベクトルと見なして、

P=A+d・t

と表せるということが基本のようです。

点 $A({x}_{A}\ \ {\ y}_{A}\ \ \ {z}_{A})$ と点 $B({x}_{B}\ \ {y}_{B}\ \ {z}_{B})$ を通る直線の方向ベクトルdは、

$d({x}_{B}-{x}_{A}\ \ {\ \ \ \ \ \ y}_{B}-{y}_{A}\ \ \ \ \ \ \ \ {\ z}_{B}-{z}_{A})$

です。

${x}_{p}={x}_{A}+{d}_{x}\bullet{}t$

${y}_{p}={y}_{A}+{d}_{y}\bullet{}t$

${z}_{p}={z}_{A}+{d}_{z}\bullet{}t$

この3式で直線が描けるわけですが、tに対して点Ｐがどんな間隔になるかは点A、B間の距離に依存します。dを単位ベクトルにすれば t=1 は長さ1に相当します。

点A、B間の距離 D は、

$D=\sqrt{{({x}_{B}-{x}_{A})}^{2}+{({y}_{B}-{y}_{A})}^{2}+{({\ z}_{B}-{z}_{A})}^{2}}\$

方向ベクトル d の単位ベクトル u は、

$u(\frac{{d}_{x}\ }{D}$ $\frac{{d}_{y}}{D}$ $\frac{{\ d}_{z}}{D}$ ) 、つまり、 $u(\frac{{x}_{B}-{x}_{A}\ }{D}$ $\frac{{y}_{B}-{y}_{A}}{D}$ $\frac{{\ z}_{B}-{z}_{A}}{D}$ )

t = 0 のとき点Pは、点 $A({x}_{A}\ \ {\ y}_{A}\ \ \ {z}_{A})$ を示し、t=D のとき点 $B({x}_{B}\ \ {y}_{B}\ \ {z}_{B})$ です。

$P({x}_{A}+{u}_{x}\bullet{}t\ \ \ {\ \ \ y}_{A}+{u}_{y}\bullet{}t\ \ \ \ \ \ {z}_{A}+{u}_{z}\bullet{}t)=$ $\ ({x}_{B}\ \ {y}_{B}\ \ {z}_{B})$

です。

Line <- function(A,B){
d <- B - A; D <- sqrt(sum(d^2))
list(d=d, D=D, u=d/D,
x=c(A[1],B[1]),y=c(A[2],B[2]),z=c(A[3],B[3]))
}

DrawLine2p <- function(A,B,aux=0,lwd=1,col="black"){
L <- Line(A,B)
lines3d(L$x,L$y,L$z,lwd=lwd,col=col)
if(aux>0){
lwd=lwd/2
A <- A - L$u*aux; B <- B + L$u*aux
lines3d(c(A[1],B[1]),c(A[2],B[2]),c(A[3],B[3]),
lwd=lwd,col=col)
}
invisible(L)
}

2点間に直線を引くことは、lines3d()を使えば良いのですが、作図では、直線がずっと続いていることを示すのに2点より長い線が適切です。

Line() は、前述の計算をします。DrawLine2p() は、auxだけ長い線を引きます。

library("rgl")
clear3d();light3d(theta=35)
spheres3d(0,0,0,radius=1,col="cyan",alpha=0.2)
DrawLine2p(c(-1,0,0),c(1,0,0),aux=0.5); text3d(1.5,0,0,"X")
DrawLine2p(c(0,-1,0),c(0,1,0),aux=0.5,col="blue"); text3d(0,1.5,0,"Y")
DrawLine2p(c(0,0,-1),c(0,0,1),aux=0.5,col="red"); text3d(0,0,1.5,"Z")
A <- c(cos(pi/6)*cos(pi/4),cos(pi/6)*sin(pi/4),sin(pi/6))
B <- c(cos(pi/3)*cos(3*pi/4),cos(pi/3)*sin(3*pi/4),sin(pi/3))
L <- DrawLine2p(A,B,lwd=2,aux=0.3);
spheres3d(L$x,L$y,L$z,radius=0.05,col=c("blue","red"),size=5)
texts3d(L$x,L$y,L$z+0.2,c("A","B"),col=c("blue","red"))
DrawLine2p(c(0,0,0),A,aux=0.3,col="green")
DrawLine2p(c(0,0,0),B,aux=0.3,col="green")
planes3d(0,0,1,col="bisque",alpha=0.9)

点が直線上にあるかどうかは、点の座標を代入して連立方程式を満たす実数 t があるかどうかで判断することになります。

2つの直線の交点は、後で考えます。平面は直線と1点で定まります。交点を持つ２直線は、同じ平面にあります。

2. 平面

平面の式を、とします。

平面は3点で決まります。3点を $({x}_{1},{y}_{1},{z}_{1}),({x}_{2},{y}_{2},{z}_{2}),({x}_{3},{y}_{3},{z}_{3})$ とします。

この点の座標を代入して、3式の連立方程式を解きます。

分母は同じでDとします。分子のOの係数を、それぞれ、 ${N}_{L}$ 、 ${N}_{M}$ 、 ${N}_{N}$ と表します。

平面の式は、

$\frac{{N}_{L}O}{D}x+\frac{{N}_{M}O}{D}y+\frac{{N}_{N}O}{D}z+O=0$

両辺にD/Oを掛けて、

${N}_{L}x+{N}_{M}x+{N}_{N}x+D=0$

したがって、3点を $({x}_{1},{y}_{1},{z}_{1}),({x}_{2},{y}_{2},{z}_{2}),({x}_{3},{y}_{3},{z}_{3})$ を含む平面の式の係数は、

${L=N}_{L}\ \ \ \ ,\ \ \ M={N}_{M}\ \ \ \ ,\ \ {\ N=N}_{N}\ \ \ ,\ \ \ O=D$

${N}_{L}$ は、

${N}_{L}=-({z}_{1}({y}_{3}-{y}_{2})-{z}_{2}{y}_{3}+{y}_{2}{z}_{3}+{y}_{1}({z}_{2}-{z}_{3}))$

$=-{y}_{3}{z}_{1}+{y}_{2}{z}_{1}+{y}_{3}{z}_{2}-{y}_{2}{z}_{3}-{y}_{1}{z}_{2}+{y}_{1}{z}_{3}$

$={y}_{3}{z}_{2}-{y}_{3}{z}_{1}-{y}_{1}{z}_{2}\ \ \ -{y}_{2}{z}_{3}+{y}_{2}{z}_{1}+{y}_{1}{z}_{3}$

一方、

h $({y}_{3}-{y}_{1})({z}_{2}-{z}_{1})-({y}_{2}-{y}_{1})({z}_{3}-{z}_{1})$

$={y}_{3}{z}_{2}-{y}_{3}{z}_{1}-{y}_{1}{z}_{2}+{y}_{1}{z}_{1}\ \ \ -{y}_{2}{z}_{3}+{y}_{2}{z}_{1}+{y}_{1}{z}_{3}-{y}_{1}{z}_{1}$

$={y}_{3}{z}_{2}-{y}_{3}{z}_{1}-{y}_{1}{z}_{2}\ \ \ -{y}_{2}{z}_{3}+{y}_{2}{z}_{1}+{y}_{1}{z}_{3}$

結局、

${L\ =N}_{L}\ \ =({y}_{3}-{y}_{1})({z}_{2}-{z}_{1})-({y}_{2}-{y}_{1})({z}_{3}-{z}_{1})$

$M={N}_{M}\ =({z}_{3}-{z}_{1})({x}_{2}-{x}_{1})-({z}_{2}-{z}_{1})({x}_{3}-{x}_{1})$

${\ N=N}_{N}\ =({x}_{3}-{x}_{1})({y}_{2}-{y}_{1})-({x}_{2}-{x}_{1})({y}_{3}-{y}_{1})$

$\ O=D\ \ \ \ \ \ {=\ -(N}_{L}\bullet{}{x}_{1}+{N}_{M}\bullet{}{y}_{1}+{N}_{N}\bullet{}{z}_{1})$

(2.1)

Rの場合、座標を c(x,y,z) のように表すのが簡便です。

p1<-c(0.4,-0.4,0.1); p2<-c(0.4,0.4,0.4); p3<-c(-0.4,-0.4,0.8)

このように座標を記述するとして平面の式の係数を求める関数は以下のようになります。

Plane <- function(p1,p2,p3){
L <- (p3[2] - p1[2]) * (p2[3] - p1[3]) - (p2[2] - p1[2]) * (p3[3] - p1[3])
M <- (p3[3] - p1[3]) * (p2[1] - p1[1]) - (p2[3] - p1[3]) * (p3[1] - p1[1])
N <- (p3[1] - p1[1]) * (p2[2] - p1[2]) - (p2[1] - p1[1]) * (p3[2] - p1[2])
O <- -(L * p1[1] + M * p1[2] + N * p1[3]);
list(L=L,M=M,N=N,O=O)
}

library("rgl")
clear3d(); Axes3d()
A<-c(0.4,-0.4,0.1); B<-c(0.4,0.4,0.4); C<-c(-0.4,-0.4,0.8)
DrawLine2p(A,B);
DrawLine2p(B,C,col="blue");
DrawLine2p(C,A,col="red");
text3d(A[1],A[2],A[3]+0.3,"A")
text3d(B[1],B[2],B[3]+0.3,"B")
text3d(C[1],C[2],C[3]+0.3,"C")
planes3d(0,0,-1,col="blue",alpha=0.7)
P <- Plane(A,B,C)
planes3d(P$L,P$M,P$N,P$O,col="cyan",alpha=0.5)
d <- sqrt(P$L^2+P$M^2+P$N^2)
a <- P$L/d; b <- P$M/d; c <- P$N/d
DrawLine2p(C,C+c(a,b,c),col="red")
DrawLine2p(C,C+c(a,0,0),col="blue")
DrawLine2p(C,C+c(0,b,0),col="blue")
DrawLine2p(C,C+c(0,0,c),col="blue")
DrawLine2p(C+c(a,b,c),C+c(a,0,0),col="green")
DrawLine2p(C+c(a,b,c),C+c(0,b,0),col="green")
DrawLine2p(C+c(a,b,c),C+c(0,0,c),col="green")

平面の式を、と表すと、係数 L、M、N、O は、平面を描画する planes3d() の引数になります。方向ベクトル(L、M、N)を点Cを起点に描くと、平面の垂線のようです。

平面は3点で決まりますが、1点と法線でも決まります。法線が同じ平面は平行な平面ですが、1点を指定すれば決まります。

この考えに従うと、法線ベクトルは平面の裏表に描けて、1つの平面の式は2つあることになります。

3点を $A({x}_{1},{y}_{1},{z}_{1}),B({x}_{2},{y}_{2},{z}_{2}),C({x}_{3},{y}_{3},{z}_{3})$ を含む平面の式の係数を点Aを通る垂線で考えてみます。

ベクトルABは、 $({x}_{2}-{x}_{1}\ \ ,\ {y}_{2}-{y}_{1}\ \ ,\ {z}_{2}-{z}_{1})$

ベクトルACは、 $({x}_{3}-{x}_{1}\ \ ,\ {y}_{3}-{y}_{1}\ \ ,\ {z}_{3}-{z}_{1})$

同様に、点B、Cを通る垂線について、同一平面上の2つのベクトルをそれぞれ表せます。

その外積を求めて、法線ベクトルを計算すれば3通りの結果になります。

点Aを基準に考えた1つ目は、連立方程式を解いたものと符号が逆になっています。両者は同じ法線を逆の向きに表しています。

2つ目と3つ目は、展開すると同じ式になります。

この結果は、連立方程式を解いたものと同じです。

3. 方向余弦

点 $A({x}_{A}\ \ {\ y}_{A}\ \ \ {z}_{A})$ と点 $B({x}_{B}\ \ {y}_{B}\ \ {z}_{B})$ の方向ベクトル (a,b,c) は2通りに計算できます。

$(a\ ,b\ ,\ c\ )=({x}_{B}-{x}_{A}\ \ {\ \ \ \ \ \ y}_{B}-{y}_{A}\ \ \ \ \ \ \ \ {\ z}_{B}-{z}_{A})$

$(a',b',c')=({x}_{A}-{x}_{B}\ \ {\ \ \ \ \ \ y}_{A}-{y}_{B}\ \ \ \ \ \ \ \ {\ z}_{A}-{z}_{B})$

数学の問題には「方向余弦を求めよ」というのがあるようです。その答えは、

$\pm{}a,\ \pm{}b,\ \pm{}c$

のように表すようです。

2つの直線の方向ベクトルが等しい場合、2直線は平行で交点がありません。2つの方向ベクトルの各要素の和が、すべてゼロになる場合も同様です。

2つの平面の法線ベクトルが等しい場合、2平面は平行で交差しません。2つの法線ベクトルの各要素の和が、すべてゼロになる場合も同様です。

4. 平面の交線

平行でない2平面は、どこかで交差します。2平面の式を以下のように表して、連立方程式を解きます。

x と、y、z の関係は以下のようになります。

したがって、2つの平面の交線は、x=0 のとき、

(4.1)

を通ります。x=1 のとき、

を通ります。2点の差は、交線の方向ベクトルです。

定数を乗じても方向に変わりはないので、分母の値を乗じると、

(4.2)

となります。

ところで、2つの平面の係数 $({L}_{1}\ \ {M}_{1}\ \ {N}_{1}),\ ({L}_{2}\ \ {M}_{2}\ \ {N}_{2})$ をベクトルと見なした外積を (cx,cy,cz) とすると、(cx,cy,cz) は以下のようになります。

(cx,cy,cz) は、方向ベクトルと同じです。

平行でない異なる2つの平面の交線は、方向ベクトルが (cx,cy,cz) で、

$(\ 0,\ \ \frac{{N}_{1}{O}_{2}-{O}_{1}{N}_{2}}{cx},\ \ \ \frac{{O}_{1}{M}_{2}-{M}_{1}{O}_{2}}{cx})$

を通ります。

ただし、2つの平面の一方がYZ平面に平行な場合は、x を変化させて y、z を求めることができません。YZ平面に平行な平面の法線は (L M N)=(L 0 0) なので、上の式の分母がゼロになります。

xc=0 の場合は、y に対して x、z を計算すれば、同様に計算できます。さらに、cy=0

であれば、z に対して x、y を計算します。

いずれの場合も、方向ベクトルは、外積を計算した (cx,cy,cz) となります。

y、z から残りを計算する連立方程式の解は、それぞれ、以下のようです。

y に対して x、z を計算を計算すると、y=0のとき交線は、

$(\frac{{O}_{1}{N}_{2}-{N}_{1}{O}_{2}}{cy},\ \ 0,\ \ \frac{{L}_{1}{O}_{2}-{O}_{1}{L}_{2}}{cy})$

を通ります。

z に対して x、y を計算を計算すると、z=0のとき交線は、

$(\frac{{{M}_{1}{O}_{2}-O}_{1}{M}_{2}}{cz},\ \ \frac{{{O}_{1}{L}_{2}-L}_{1}{O}_{2}}{cz}\ ,\ \ 0)$

を通ります。

5. 平面の成す角

平行でない異なる2つの平面は、どこかで交差し、交差角は法線同士の交差角度です。

2つの平面を以下とすると、

2つの平面のなす角 θ は、法線の方向ベクトルを、それぞれベクトル a、b として、以下のように計算されます。


	(5.1)

これは、ベクトルa、b の単位ベクトルの内積を計算しています。ベクトルa、b が元々単位ベクトルなら分子が cosθ です。したがって、値域は [-1,1] です。

また、異なる3点から求めた平面の式の係数、Ｌ、Ｍ、Ｎはすべてがゼロにはならないので、分母がゼロになることはありません。

また、cosθ の値域は、[-1,1] で、その逆関数は[0,π] を返します。cos(π/4)、cos(-π/4)、cos( $\frac{7π}{4}$ ) はすべて $\frac{\sqrt{2}}{2}$ なので区別はありません。

2つの平面のなす角は、4か所測ることができ、直交している場合を除き、2つの値になります。2つの値の和はπです。空間に目印がなければ、角度は正の値ではかられ、小さい方を取るなどと決めることになります。この場合、平面のなす角の値域は [0,π/2] です。

6. 卯酉線の例

XY平面を赤道面、Z軸を自転軸とした地球楕円体を考えます。ただし、赤道径を1、極半径 b を b=0.6 とした扁平な地球を描きます。

緯度経度が (φ、λ) の地点Aを通る垂線を引き、その垂線を含んで、λの経線が作る面に直交する平面を描きます。赤道面と、この平面がなす角度は φ です。この平面による楕円体の断面の周は卯酉線です。 (φ、λ) の地点Aにおいて経線（南北）に直交する線で、真東、真西を延長したものです。

「１．地球楕円体上の点」に従って、 A(φ、λ) の座標 A $({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})$ は以下のようになります。

$\beta{}=atan(b\cdot{}tan\varphi{})$

$({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})=(cos\beta{}\cdot{}cos\lambda{}\ ,\ \ cos\beta{}\cdot{}sin\lambda{}\ ,\ \ b\cdot{}sin\beta{})$ (6.1)

経度 λ の経線の作る平面を新たにX'Z平面とします。この平面でAの座標は、

A(cosβ、b・sinβ)

です。Aで地面に垂直な直線AQは、赤道面と φ で交差します。この直線の切片を q とすると、

$y=tan\varphi{}\cdot{}x+q$

と表され、点Aの座標を代入して切片は、

$q=b\cdot{}sin\beta{}-cos\beta{}\cdot{}tan\varphi{}$ (6.2)

となります。点Sの座標は、S(0, q) です。

経度 λ の経線の作る平面は、点Aを含み、原点、両極を含みます。

$({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})$ 、(0,0,0)、(0,0,1) を含む平面は、

${y}_{A}\cdot{}x-{x}_{A}\cdot{}y=0$ (6.3)

で、これに直交する平面の法線ベクトル (a b c) は ${y}_{A}a-{x}_{A}\bullet{}b+0\bullet{}c=0$ を満たすので、 $(a\ b\ c)=({x}_{A}\ \ \ {y}_{A}\ \ \ c)$ として、経度 λ の経線の作る平面に直交する平面は、

${x}_{A}\cdot{}x+{y}_{A}\cdot{}y+c\cdot{}z+O=0$

(6.4)

この平面は、点A、Sを含むので、 $({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})$ 、(0 0 q) を代入して、

したがって、卯酉線の描く平面は、

(6.5)

緯度経度(φ=π/6、λ=π/4) について描いて見ます。

b <- 0.6
φ <- pi/6; λ <- pi/4
DrawEarthEllipsoid(b=b,drawspheroid=FALSE,alpha=0.7)
β <- atan(b*tan(φ))
xp <- cos(β)*cos(λ); yp <- cos(β)*sin(λ); zp <- b*sin(β)
spheres3d(xp,yp,zp,radius=0.05,col="red")
q <- b*sin(β)-cos(β)*tan(φ)
DrawLine2p(c(0,0,q),c(xp,yp,zp))
O <- -((yp^2+xp^2)*q)/(q-zp)
c <- (yp^2+xp^2)/(q-zp)
planes3d(yp,-xp,0,0,col="magenta",alpha=0.6)
planes3d(xp,yp,c,O,col="green",alpha=1)
DrawLine2p(c(-cos(λ),-sin(λ),0),c(cos(λ),sin(λ),0),lwd=3)

7. 赤道面となす角

卯酉線を描く平面 (xp,py,c,O) と、赤道面 (0,0,1) のなす角 θ は、

$cos\theta{}=\frac{c}{\sqrt{{xp}^{2}+{yp}^{2}+{c}^{2}}\ }$ (6.6)

赤道面を (0,0,-1)と表せば、

$cos\theta{}=\frac{ーｃ}{\sqrt{{xp}^{2}+{yp}^{2}+{c}^{2}}\ }$

です。θは、π/6 か 5π/6 のいずれかになります。

平面の式は、3点を与える順番で変わります。

また、南緯のつもりで φ=-π/6 を計算しても、π/6 か 5π/6 のいずれかになることに変わりはありません。

平面のなす角は、符号や90°より大きい小さいといったことで、方向を判断することには役立たないようです。

8. 赤道面との交線

前述の、

$({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})=(cos\beta{}\cdot{}cos\lambda{}\ ,\ \ cos\beta{}\cdot{}sin\lambda{}\ ,\ \ b\cdot{}sin\beta{})$ (6.1)

$q=b\cdot{}sin\beta{}-cos\beta{}\cdot{}tan\varphi{}$ (6.2)

から、鉛直線の式は、

$x={x}_{A}+{(0-x}_{A})\cdot{}t$

$y={y}_{A}+{(0-y}_{A})\cdot{}t$

$z={z}_{A}+{(q-z}_{A})\cdot{}t$

と表せます。 $({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})\$ から (0 0 q) へ向かう方向ベクトルは、 $({-x}_{A}\ \ -{y}_{A}\ \ \ q-{z}_{A})\$ です。

A点における鉛直線と赤道面との交点Qは、z = 0 となる t を求めると、

${z}_{A}+{(q-z}_{A})\cdot{}t=0$

$t=-{z}_{A}\ /\ (q-{z}_{A})$

となり、

Q $({x}_{Q}\ \ {\ y}_{Q}\ \ 0)=(q\cdot{}{x}_{A}/{(q-z}_{A})\ \ \ \ \ \ \ \ q\cdot{}{y}_{A}/{(q-z}_{A})\ \ \ \ \ \ 0)$ (6.2.1)

卯酉線を描く平面の式(6.5) は、z=0 のとき、

(10.6.2)

となって、卯酉線を描く平面と、赤道面の交線となります。

x軸との交点は、

$({x}_{c}\ ,\ 0,\ 0)=(\frac{({{x}_{A}}^{2}+{{y}_{A}}^{2})q}{(q-{z}_{A}){x}_{A}}\ \ ,\ \ 0,\ \ 0)$

y軸との交点は

$(0,\ {y}_{c},\ 0)=(0\ ,\ \ \frac{({{x}_{A}}^{2}+{{y}_{A}}^{2})q}{(q-{z}_{A}){y}_{A}}\ \ ,\ \ 0)$

卯酉線と赤道の交点は、式(10.6.2) と、単位円の交点です。あるいは、卯酉線を描く平面と、単位円の交点です。

$u={{x}_{A}}^{2}+{{y}_{A}}^{2}\ \ ,\ \ \ v=q-{z}_{A}\$ と置いて、

9. 楕円体の截線

「断面」は「面」と聞こえるので、その周は「切り（截り）線」と言うようです。

卯酉線は垂直截線で、地球楕円体の表面と平面の交線です。楕円体と、平面の式は、

地点Aの緯度経度(φ、λ)は、赤道面をXY平面としたXYZ直交座標で、

$\beta{}=atan(b\bullet{}tan\varphi{})$

$({x}_{A}\ \ \ {y}_{A}\ \ \ {z}_{A})=(cos\beta{}\bullet{}cos\lambda{},\ cos\beta{}\bullet{}sin\lambda{},\ b\bullet{}sin\beta{})$

長半径が1、短半径がb の地球楕円体で、地点A(φ、λ)において、鉛直線は(0,0,q)を通ります。

$q=b\bullet{}sin\beta{}-cos\beta{}\bullet{}tan\varphi{}$

$u=\frac{{{x}_{A}}^{2}+{{y}_{A}}^{2}}{q-{z}_{A}}$ と置いて、平面の式から z を求めて、楕円体の式に代入すると、

${x}^{2}+{y}^{2}+\frac{{({x}_{A}\cdot{}x+{y}_{A}\cdot{}y-qu)}^{2}}{{b}^{2}{u}^{2}}=1$

${b}^{2}{u}^{2}$ を両辺に乗じて整理すると、

(6.3.1)

これは楕円を含む2次曲線（円錐曲線）の一般形です。楕円体の平面による截線のXY平面への投影です。一般形は、

$A{x}^{2}+Bxy+C{y}^{2}+Dx+Ey+F=0$

(6.3.2)

と表され、その係数A、B、C、D、E、F から楕円の標準形を導いて、一周分のｘ、ｙ座標を算出し、それを平面の式に代入して、zを計算すれば卯酉線が引けます。

「４．卯酉線」、「５．楕円の標準形」

b <- 0.6 # 赤道半径1、極半径bの準拠楕円体
φ <- pi/6; λ <- pi/4 # 地点A 緯度30°、経度45°
β <- atan(b*tan(φ)) # 地点Aの更成緯度
xA <- cos(β)*cos(λ); yA <- cos(β)*sin(λ)
zA <- b*sin(β)
q <- b*sin(β)-cos(β)*tan(φ) # 垂線の切片
DrawEarthEllipsoid(b=b,drawspheroid=FALSE,alpha=0.2)
spheres3d(xA,yA,zA,radius=0.05,col="red") # 地点A
DrawLine2p(c(0,0,q),c(xA,yA,zA)) # 垂線
LongitudeLine(λ,b) # 経度λの経線
# 卯酉線の赤道面への投影図形の式の係数
u <- (yA^2+xA^2)/(q-zA)
A <- xA^2+b^2*u^2; B <- 2*xA*yA
C <- yA^2+b^2*u^2
D <- -2*xA*q*u; E <- -2*yA*q*u
F <- u^2*(q^2-b^2)
u <- StandardForm(A,B,C,D,E,F) # 標準形の諸元
t <- seq(0,2*pi,length.out=361)
x <- u$a*cos(t); y <- u$b*sin(t) # 標準形の周
g <- u$vectors %*% rbind(x,y) # 回転
x2 <- g[1,]+u$cx; y2 <- g[2,]+u$cy # 平行移動
lines3d(x2,y2,rep(0,length(x)),col="orange")
z <- -(xA*x2+yA*y2-R*q)/R # x,yからz算出
lines3d(x2,y2,z,lwd=1,col="blue") # 截線描画
for(i in seq(1,361,by=15))
DrawLine2p(c(x2[i],y2[i],0),c(x2[i],y2[i],z[i]))