adsomath.h

Go to the documentation of this file.

/*****************************************************************************
 * $Id: adsomath.h,v 1.16 2003/06/27 08:36:41 schou Exp $
 * Author: Jakob Schou Jensen <schou@sourceforge.net>
 * Licence: GPL
 *
 * Description: Adso math
 *
 ****************************************************************************/


#ifndef ADSOMATH_H
#define ADSOMATH_H

#include <stdlib.h>
#include <math.h>
#include <adsobase.h>

#ifdef ADSO_WINDOWS
#define M_PI    3.14159265358979323846
#define M_TWOPI         (M_PI * 2.0)
#define M_PI_2    1.57079632679489661923
#define M_PI_4    0.78539816339744830962
#define M_3PI_4   2.3561944901923448370E0
inline float sin(float v){ return (float)sin((double)v); }
inline float cos(float v){ return (float)cos((double)v); }
inline float tan(float v){ return (float)tan((double)v); }
inline float asin(float v){ return (float)asin((double)v); }
inline float acos(float v){ return (float)acos((double)v); }
inline float atan(float v){ return (float)atan((double)v); }
inline float sqrt(float v){ return (float)sqrt((double)v); }
inline float floor(float v){ return (float)floor((double)v); }
inline float ceil(float v){ return (float)ceil((double)v); }
inline float fabs(float v){ return (float)fabs((double)v); }
inline float hypot(float a1, float a2){ return (float)hypot((double)a1,(double)a2); }
inline float pow(float a1, float a2){ return (float)pow((double)a1,(double)a2); }
inline int isnan(float v){ return 0; } // XXX fix this
inline int isinf(float v){ return 0; } // XXX fix this
#endif

// Basic math
inline double sqr(double v){return v*v;}
inline double cub(double v){return v*v*v;}
inline float sqr(float v){return v*v;}
inline float cub(float v){return v*v*v;}

// Clamp to [-1,1]
inline float clamp(float v){
  if(v<-1.0) return -1.0;
  else if(v>1.0) return 1.0;
  return v;
}

inline double clamp(double v){
  if(v<-1.0) return -1.0;
  else if(v>1.0) return 1.0;
  return v;
}

// PI normalize [-PI,PI]
inline float pin(float v){
  while(v>M_PI) v-=(float)M_PI*2.0f;
  while(v<-M_PI) v+=(float)M_PI*2.0f;
  return v;
}

inline double pin(double v){
  while(v>M_PI) v-=M_PI*2.0;
  while(v<-M_PI) v+=M_PI*2.0;
  return v;
}

inline float sign(float v){
  if(v<0.0) return -1.0;
  else return 1.0;
}

inline double sign(double v){
  if(v<0.0) return -1.0;
  else return 1.0;
}

// frandom(mag) : returns a random value in the intervaæ ]mag:0]
#ifdef ADSO_WINDOWS
//inline unsigned int random(){ return (unsigned int)rand() * (unsigned int)rand(); }
//inline void srandom(unsigned int seed){ srand(seed); }
unsigned int random();
void srandom(unsigned int seed);
double frandom(double mag);
inline float frandom(float mag){ return (float)frandom((double)mag); }
#else
inline double frandom(double mag){
  return (double)random()/(double)RAND_MAX*mag;
}
#endif


// Vector 2D double
class v2d {
public:
  double x,y;

  v2d hat(){
    return v2d(-y,x);
  }

  double length(){
    return hypot(x,y);
  }

  double dist(v2d p){
    return hypot(x-p.x,y-p.y);
  }

  v2d norm(){
    double l=hypot(x,y);
    return v2d(x/l,y/l);
  }

  v2d operator=(const v2d v){
    x = v.x;
    y = v.y;
    return v2d(x,y);
  }

  int operator==(const v2d a){
    return x==a.x && y==a.y;
  }

  v2d operator+(const v2d a){
    return v2d(x+a.x,y+a.y);
  }

  v2d operator-(const v2d a){
    return v2d(x-a.x,y-a.y);
  }

  v2d operator*(const v2d a){
    return v2d(x*a.x,y*a.y);
  }

  v2d operator/(const double d){
    return v2d(x/d,y/d);
  }

  v2d operator*(const double m){
    return v2d(x*m,y*m);
  }

  friend v2d operator*(const double m,v2d v){
    return v2d(v.x*m,v.y*m);
  }

  v2d(double cx, double cy){
    x=cx;
    y=cy;
  }

  v2d(){
    x=0; y=0;
  }
};

inline double det(v2d a, v2d b){
  return a.x*b.y-a.y*b.x;
}

// Vector 2D float
class v2f {
public:
  float x,y;

  v2f hat(){
    return v2f(-y,x);
  }

  float length(){
    return (float)hypot(x,y);
  }

  float dist(v2f p){
    return (float)hypot(x-p.x,y-p.y);
  }

  v2f norm(){
    float l = (float)hypot(x,y);
    return v2f(x/l,y/l);
  }

  v2f operator=(const v2f v){
    x = v.x;
    y = v.y;
    return v2f(x,y);
  }

  int operator==(const v2f a){
    return x==a.x && y==a.y;
  }

  v2f operator+(const v2f a){
    return v2f(x+a.x,y+a.y);
  }

  v2f operator-(const v2f a){
    return v2f(x-a.x,y-a.y);
  }

  v2f operator*(const v2f a){
    return v2f(x*a.x,y*a.y);
  }

  v2f operator/(const float d){
    return v2f(x/d,y/d);
  }

  v2f operator*(const float m){
    return v2f(x*m,y*m);
  }

  friend v2f operator*(const float m,v2f v){
    return v2f(v.x*m,v.y*m);
  }

  v2f(float cx, float cy){
    x=cx;
    y=cy;
  }

  v2f::v2f(class p3f p);

  v2f(){
    x=0; y=0;
  }
};

inline float det(v2f a, v2f b){
  return a.x*b.y-a.y*b.x;
}

// Vector 3D float
class v3f {
public:
  float x,y,z;

  float length(){
    return (float)sqrt(x*x+y*y+z*z);
  }

  float dist(v3f p){
    return (float)sqrt(((x-p.x)*(x-p.x))+((y-p.y)*(y-p.y))+((z-p.z)*(z-p.z)));
  }

  v3f norm(){
    float l = sqrt(x*x+y*y+z*z);
    //if(l<1e-20&&0) return v3f(1.0,0.0,0.0);
    return v3f(x/l,y/l,z/l);
  }

  v3f operator=(const v3f v){
    x = v.x;
    y = v.y;
    z = v.z;
    return v3f(x,y,z);
  }

  int operator==(const v3f a){
    return x==a.x && y==a.y && z==a.z;
  }

  v3f operator+(const v3f a){
    return v3f(x+a.x,y+a.y,z+a.z);
  }

  v3f operator-(const v3f a){
    return v3f(x-a.x,y-a.y,z-a.z);
  }

  v3f operator/(const float d){
    return v3f(x/d,y/d,z/d);
  }

  float operator*(const v3f v){
    return x*v.x + y*v.y + z*v.z;
  }

  v3f operator*(const float m){
    return v3f(x*m,y*m,z*m);
  }

  friend v3f operator*(const float m,const v3f v){
    return v3f(v.x*m,v.y*m,v.z*m);
  }

  v3f operator-(){
    return v3f(-x,-y,-z);
  }

  v3f operator-=(const v3f a){
    x = x-a.x;
    y = y-a.y;
    z = z-a.z;
    return v3f(x,y,z);
  }

  v3f operator+=(const v3f a){
    x = x+a.x;
    y = y+a.y;
    z = z+a.z;
    return v3f(x,y,z);
  }

  v3f(float cx, float cy, float cz){
    x=cx;
    y=cy;
    z=cz;
  }

  v3f(class p3f p);
  v3f(class s2f p);
  v3f(class s2d p);

  v3f(){
    x=0; y=0; z=0;
  }
};

inline v3f cross(v3f a, v3f b){
  return v3f(a.y*b.z-a.z*b.y, 
       a.z*b.x-a.x*b.z, 
       a.x*b.y-a.y*b.x);
}

// Polar vector 3D float
class p3f {
public:
  float t,p,r; // theta, phi, radius (lenght)

  p3f operator*(const float m){
    return p3f(t,p,r*m);
  }

  float sphereDist(p3f q){
    return (float)(acos(cos(t)*cos(q.t)+sin(t)*sin(q.t)*cos(q.p-p)));
  }

  p3f(float ct, float cp, float cr){
    p=cp;
    t=ct;
    r=cr;
  }

  p3f(){
    p=0; t=0; r=0;
  }
};


inline v3f::v3f(p3f p){
  x = (float)(sin(p.t)*cos(p.p)*p.r);
  y = (float)(sin(p.t)*sin(p.p)*p.r);
  z = (float)(cos(p.t)*p.r);
}

// Sphere vector Float
class s2f {
public:
  float t,p; // theta, phi, radius (lenght)

  s2f(float ct, float cp){
    p=cp;
    t=ct;
  }

  s2f(){
    p=0; t=0;
  }

  s2f(const p3f v){
    t = v.t;
    p = v.p;
  }

  s2f(const v3f v){
    v3f n = v;
    float s;

    n = n.norm();
    s = (float)sqrt(1.0-sqr(n.z));

    t = (float)acos(n.z);
    if(s<1e-10) p = 0;
    else p = (float)acos(clamp(n.x/s))*sign(n.y);
  }
};

// Sphere vector Double
class s2d {
public:
  double t,p; // theta, phi, radius (lenght)

  s2d zy_rotate(double at, double ap){
    double cz,sac_cz,cx;

    cz = clamp((float)(cos(at)*cos(t)-sin(at)*sin(t)*cos(ap+p)));
    sac_cz = sqrt(1.0-sqr(cz));
    cx = sin(at)*cos(t)+cos(at)*sin(t)*cos(ap+p);

    if(sac_cz>1e-10)
      return s2d(acos(cz), acos(clamp(cx/sac_cz))*sign(pin(ap+p)));
    else
      return s2d(acos(cz), 0.0);
  }

  s2d yz_rotate(double at, double ap){
    double cz,sac_cz,cx;

    cz = clamp(cos(at)*cos(t)-sin(at)*sin(t)*cos(p));
    sac_cz = sqrt(1.0-sqr(cz));
    cx = sin(at)*cos(t)+cos(at)*sin(t)*cos(p);

    if(sac_cz>1e-10)
      return s2d(acos(cz), pin(acos(clamp(cx/sac_cz))*sign(pin(p))+ap));
    else
      return s2d(acos(cz), 0.0);
  }

  double sphereDist(s2d q){
    return acos(cos(t)*cos(q.t)+sin(t)*sin(q.t)*cos(q.p-p));
  }

  double sphereDist_dt(s2d q){
    double d = sqrt(1.0-sqr(cos(t)*cos(q.t)+sin(t)*sin(q.t)*cos(q.p-p)));
    if(d<1e-10) return 0.0;
    else return (cos(t)*sin(q.t)-sin(t)*cos(q.t)*cos(p-q.p))/d;
  }

  double sphereDist_dp(s2d q){
    double d = sqrt(1.0-sqr(cos(t)*cos(q.t)+sin(t)*sin(q.t)*cos(q.p-p)));
    if(d<1e-10) return 0.0;
    else return (sin(t)*sin(q.t)*sin(q.p-p))/d;
  }

  s2d(double ct, double cp){
    p=cp;
    t=ct;
  }

  s2d(){
    p=0; t=0;
  }
};

inline v3f::v3f(s2f p){
  x = (float)(sin(p.t)*cos(p.p));
  y = (float)(sin(p.t)*sin(p.p));
  z = (float)cos(p.t);
}

inline v3f::v3f(s2d p){
  x = (float)(sin(p.t)*cos(p.p));
  y = (float)(sin(p.t)*sin(p.p));
  z = (float)cos(p.t);
}

// used by some pseudo realistic texture coord generator
inline v2f::v2f(p3f p){
  x = (float)(sin(p.t)*cos(p.p));
  y = (float)(sin(p.t)*sin(p.p));
}

// Square Matrix 3x3 Float
class sqm3f {
 public:
  float m[3][3];

  v3f col(int i){
    return v3f(m[0][i],m[1][i],m[2][i]);
  }

  v3f row(int i){
    return v3f(m[i][0],m[i][1],m[i][2]);
  }

  void setcol(int i, v3f c){
    m[0][i] = c.x; m[1][i] = c.y; m[2][i] = c.z;
  }

  void setrow(int i, v3f c){
    m[i][0] = c.x; m[i][1] = c.y; m[i][2] = c.z;
  }

  sqm3f operator=(const sqm3f v){
    m[0][0] = v.m[0][0];  m[0][1] = v.m[0][1];  m[0][2] = v.m[0][2];
    m[1][0] = v.m[1][0];  m[1][1] = v.m[1][1];  m[1][2] = v.m[1][2];
    m[2][0] = v.m[2][0];  m[2][1] = v.m[2][1];  m[2][2] = v.m[2][2];
    return sqm3f(m);
  }

  sqm3f operator+=(const sqm3f a){
    m[0][0] += a.m[0][0];  m[0][1] += a.m[0][1];  m[0][2] += a.m[0][2];
    m[1][0] += a.m[1][0];  m[1][1] += a.m[1][1];  m[1][2] += a.m[1][2];
    m[2][0] += a.m[2][0];  m[2][1] += a.m[2][1];  m[2][2] += a.m[2][2];
    return sqm3f(m);
  }

  sqm3f operator+(const sqm3f a){
    sqm3f r;
    r.m[0][0] = m[0][0] + a.m[0][0];  r.m[0][1] = m[0][1] + a.m[0][1];  r.m[0][2] = m[0][2] + a.m[0][2];
    r.m[1][0] = m[1][0] + a.m[1][0];  r.m[1][1] = m[1][1] + a.m[1][1];  r.m[1][2] = m[1][2] + a.m[1][2];
    r.m[2][0] = m[2][0] + a.m[2][0];  r.m[2][1] = m[2][1] + a.m[2][1];  r.m[2][2] = m[2][2] + a.m[2][2];
    return r;
  }

  sqm3f operator-=(const sqm3f a){
    m[0][0] -= a.m[0][0];  m[0][1] -= a.m[0][1];  m[0][2] -= a.m[0][2];
    m[1][0] -= a.m[1][0];  m[1][1] -= a.m[1][1];  m[1][2] -= a.m[1][2];
    m[2][0] -= a.m[2][0];  m[2][1] -= a.m[2][1];  m[2][2] -= a.m[2][2];
    return sqm3f(m);
  }

  sqm3f operator-(const sqm3f a){
    sqm3f r;
    r.m[0][0] = m[0][0] - a.m[0][0];  r.m[0][1] = m[0][1] - a.m[0][1];  r.m[0][2] = m[0][2] - a.m[0][2];
    r.m[1][0] = m[1][0] - a.m[1][0];  r.m[1][1] = m[1][1] - a.m[1][1];  r.m[1][2] = m[1][2] - a.m[1][2];
    r.m[2][0] = m[2][0] - a.m[2][0];  r.m[2][1] = m[2][1] - a.m[2][1];  r.m[2][2] = m[2][2] - a.m[2][2];
    return r;
  }

  sqm3f operator*(const sqm3f a){
    sqm3f r;
    int i,j;
    for(j=0;j<3;j++){
      for(i=0;i<3;i++){
  r.m[i][j] = m[i][0]*a.m[0][j] + m[i][1]*a.m[1][j] + m[i][2]*a.m[2][j];
      }
    }
    return r;
  }

  int operator==(const sqm3f a){
    int i,j;
    for(j=0;j<3;j++){
      for(i=0;i<3;i++){
  if(a.m[i][j] != m[i][j]) return 0;
      }
    }
    return 1;
  }

  v3f operator*(const v3f v){
    return v3f(m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z,
         m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z,
         m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z);
  }

  sqm3f operator*(const float f){
    sqm3f p;
    p.m[0][0] = m[0][0] * f;  p.m[0][1] = m[0][1] * f;  p.m[0][2] = m[0][2] * f;
    p.m[1][0] = m[1][0] * f;  p.m[1][1] = m[1][1] * f;  p.m[1][2] = m[1][2] * f;
    p.m[2][0] = m[2][0] * f;  p.m[2][1] = m[2][1] * f;  p.m[2][2] = m[2][2] * f;
    return p;
  }

  sqm3f operator*=(const float f){
    m[0][0] *= f;  m[0][1] *= f;  m[0][2] *= f;
    m[1][0] *= f;  m[1][1] *= f;  m[1][2] *= f;
    m[2][0] *= f;  m[2][1] *= f;  m[2][2] *= f;
    return sqm3f(m);
  }

  sqm3f trans(){
    sqm3f r;
    r.m[0][0] = m[0][0];  r.m[0][1] = m[1][0];  r.m[0][2] = m[2][0];
    r.m[1][0] = m[0][1];  r.m[1][1] = m[1][1];  r.m[1][2] = m[2][1];
    r.m[2][0] = m[0][2];  r.m[2][1] = m[1][2];  r.m[2][2] = m[2][2];    
    return r;
  }

  sqm3f inverse();

  // Vector square: a*b^T
  void vsqr(v3f a, v3f b){
    m[0][0] = a.x*b.x;  m[0][1] = a.x*b.y;  m[0][2] = a.x*b.z;
    m[1][0] = a.y*b.x;  m[1][1] = a.y*b.y;  m[1][2] = a.y*b.z;
    m[2][0] = a.z*b.x;  m[2][1] = a.z*b.y;  m[2][2] = a.z*b.z;
  }

  void ident(){
    m[0][0] = 1.0;  m[0][1] = 0.0;  m[0][2] = 0.0;
    m[1][0] = 0.0;  m[1][1] = 1.0;  m[1][2] = 0.0;
    m[2][0] = 0.0;  m[2][1] = 0.0;  m[2][2] = 1.0;
  }

  void zero(){
    m[0][0] = 0.0;  m[0][1] = 0.0;  m[0][2] = 0.0;
    m[1][0] = 0.0;  m[1][1] = 0.0;  m[1][2] = 0.0;
    m[2][0] = 0.0;  m[2][1] = 0.0;  m[2][2] = 0.0;
  }

  sqm3f(float a[3][3]){
    m[0][0] = a[0][0];  m[0][1] = a[0][1];  m[0][2] = a[0][2];
    m[1][0] = a[1][0];  m[1][1] = a[1][1];  m[1][2] = a[1][2];
    m[2][0] = a[2][0];  m[2][1] = a[2][1];  m[2][2] = a[2][2];
  }

  sqm3f(v3f a, v3f b, v3f c){
    m[0][0] = a.x;  m[0][1] = b.x;  m[0][2] = c.x;
    m[1][0] = a.y;  m[1][1] = b.y;  m[1][2] = c.y;
    m[2][0] = a.z;  m[2][1] = b.z;  m[2][2] = c.z;
  }

  sqm3f(){
    m[0][0] = 0.0;  m[0][1] = 0.0;  m[0][2] = 0.0;
    m[1][0] = 0.0;  m[1][1] = 0.0;  m[1][2] = 0.0;
    m[2][0] = 0.0;  m[2][1] = 0.0;  m[2][2] = 0.0;
  }
};

inline sqm3f rotateMatrix(float a, v3f v){
  sqm3f U,S,R;

  U.vsqr(v,v);
  S.m[0][0] = 0.0; S.m[0][1] = -v.z; S.m[0][2] = v.y;
  S.m[1][0] = v.z; S.m[1][1] = 0.0; S.m[1][2] = -v.x;
  S.m[2][0] = -v.y; S.m[2][1] = v.x; S.m[2][2] = 0.0;

  R.ident();
  R -= U;
  R *= (float)cos(a);
  R += U;
  S *= (float)sin(a);
  R += S;
  
  return R;
}

// rotate u a around v
inline v3f rotate(float a, v3f v, v3f u){
  sqm3f U,S,R;

  U.vsqr(v,v);
  S.m[0][0] = 0.0;
  S.m[0][1] = -v.z;
  S.m[0][2] = v.y;
  S.m[1][0] = v.z;
  S.m[1][1] = 0.0;
  S.m[1][2] = -v.x;
  S.m[2][0] = -v.y;
  S.m[2][1] = v.x;
  S.m[2][2] = 0.0;

  R.ident();
  R -= U;
  R *= cos(a);
  R += U;
  S *= sin(a);
  R += S;
  
  return R*u;
}

inline sqm3f star(v3f v){
  sqm3f r;
  r.m[0][0] = 0;    r.m[0][1] = -v.z; r.m[0][2] = v.y;
  r.m[1][0] = v.z;  r.m[1][1] = 0;    r.m[1][2] = -v.x;
  r.m[2][0] = -v.y; r.m[2][1] = v.x;  r.m[2][2] = 0;    
  return r;
}



class line2d {
public:
  v2d v0, v1;

  int intersect(line2d l);
  line2d(){}
  line2d(v2d p0, v2d p1){
    v0 = p0;
    v1 = p1;
  }
};

class tri2d {
public:
  v2d v0, v1, v2;

  tri2d(){}
  tri2d(v2d p0, v2d p1, v2d p2){
    v0 = p0;
    v1 = p1;
    v2 = p2;
  }
};

// 3D triangle
class tri3f {
public:
  v3f v0, v1, v2;

  tri3f(){}
  tri3f(v3f a, v3f b, v3f c){ v0=a; v1=b; v2=c;}
};

// 3D triangle with surface normals
class tri3n {
public:
  v3f v0, v1, v2;
  v3f n0, n1, n2;

  tri3n(){}
};

// 3D vector and surface normal
class vn3f {
public:
  v3f v;
  v3f n;

  vn3f(){}
  vn3f(v3f sv, v3f sn){v=sv;n=sn;}
};

// 3D vector, surface normal, and texture coordinate
class vnt3f {
public:
  v3f v;
  v3f n;
  v2f t;

  vnt3f(){}
  vnt3f(v3f sv, v3f sn, v2f st){v=sv;n=sn;t=st;}
};

// 3D vector, surface normal, and two texture coordinate sets
class vntt3f {
public:
  v3f v;
  v3f n;
  v2f t1,t2;

  vntt3f(){}
  vntt3f(v3f sv, v3f sn, v2f st1, v2f st2){ v=sv; n=sn; t1=st1; t2=st2; }
};


/* simple set (may contain doubles) */
template<class T>
class Set {
  int array_size;
  int count;
  T *array;
public:
  T get(int i){ return array[i]; }

  /*
  void add(tri_set *ts){
    int i,s = tris+ts->count();
    if(s>array_size){
      if(array_size*2>s) s = array_size*2;
      tri *n = new tri[s];
      for(i=0;i<tris;i++){
  n[i] = array[i];
      }
      delete(array);
      array = n;
      array_size = s;
    }
    for(i=0;i<ts->count();i++)
      array[tris+i] = ts->get(i);
    tris += ts->count();
  }
  */

  void add(T t){
    if(count==array_size){
      T *n = new T[array_size*2];
      int i;
      for(i=0;i<array_size;i++){
  n[i] = array[i];
      }
      delete(array);
      array = n;
      array_size *= 2;
    }
    array[count++] = t;
  }

  int size(){ return count; }

  void clear(){
    delete(array);
    array=new T[16];
    array_size=16;
    count=0;
  }

  Set(){
    array=new T[16];
    array_size=16;
    count=0;
  }
};

class SphereLine {
  p3f a,b,c;

 public:
  p3f linePoint(float t){
    float x,z,cp;
    z = clamp(cos(a.t)*cos(t*c.t)-sin(a.t)*sin(t*c.t)*cos(c.p));
    x = sin(a.t)*cos(t*c.t)+cos(a.t)*sin(t*c.t)*cos(c.p);
    if(1.0-z<1e-10) cp = 0;
    else cp = clamp(x/sqrt(1.0f-sqr(z)));
    
    return p3f( acos(z), pin(acos(cp)*sign(c.p)+a.p), 1.0 );
  }

  void generateLine(Set<p3f> *ps, int count){
    int j;
    for(j=0;j<count;j++){
      p3f lp = linePoint((float)j/(float)(count-1));
      ps->add(lp);
    }
  }

  void generateLine(p3f *pl, int count){
    int j;
    for(j=0;j<count;j++){
      pl[j] = linePoint((float)j/(float)(count-1));
    }
  }

  void endPoints(p3f la, p3f lb){
    float z,x,cp;

    a = la;
    b = lb;
    
    z = clamp(cos(a.t)*cos(b.t) + sin(a.t)*sin(b.t)*cos(a.p-b.p));
    c.t = acos(z);

    x = cos(a.t)*sin(b.t)*cos(a.p-b.p) - sin(a.t)*cos(b.t);
    if(1.0-z<1e-10) cp = 0.0;
    else cp = x/sqrt(1.0f-sqr(z));
    c.p = acos(cp) * sign(pin(b.p-a.p));
  }

  SphereLine(){}
};


#endif

---