awk一行代码将平面在极折赤平投影上

awk一行代码在极折赤平投影上画平面

最近用gmt6在极折赤平投影上投点和线的时候遇到了一些问题，于是把脚本贴上来。
输入的三列数据为平面的法向量坐标

gmt begin test png
gmt basemap -JX2c/2c -R-10/10/-10/10 -Baf 
awk 'BEGIN{
    pid=0.017453
    for (az=0; az<=360; az++) {
    print 10*cos(az*pid),10*sin(az*pid) }
}' | gmt clip -W1.0p,black 
echo 0 1 1 | awk '{
    xn=$1;          yn=$2;          zn=$3;
    x0=0;           y0=0;           z0=1;
    xp=yn*z0-y0*zn; yp=x0*zn-xn*z0; zp=xn*y0-x0*yn;
    xq=yn*zp-yp*zn; yq=xp*zn-xn*zp; zq=xn*yp-xp*yn;
    print ">"
    for (i=0; i<360; i+=5) {
        xl=xp*cos(i*0.017453)+xq*sin(i*0.017453)
        yl=yp*cos(i*0.017453)+yq*sin(i*0.017453)
        zl=zp*cos(i*0.017453)+zq*sin(i*0.017453)
        azi=atan2(xl,yl)
        if (zl<0) azi=3.1416+azi
        thi=atan2(sqrt(zl*zl),sqrt(xl*xl+yl*yl))
        x=10*cos(thi)*sin(azi)
        y=10*cos(thi)*cos(azi)
        if (zl>0) print x,y
    }
}' | gmt plot -W0.3p,red
gmt clip -C
gmt end

原理是平面的法向量n(xn,yn,zn)先和001向量叉乘得到一个该平面的一个基向量p(xp,yp,zp)，然后该基向量p再和法向量n叉乘得到该平面的另一个基向量q(xq,yq,zq),然后以pq两个基向量为基底画圆即可。
效果如下图