#if !defined(__VECTORN_H)
#define __VECTORN_H
#include  <stdio.h>
//プリコンパイルヘッダが作成できないメッセージを抑制する
#pragma option -w-pch
//-------------------
/*
  ★警告　すべてのモジュールとその相互関係が試験できているわけではない。
  ★c:\src\commmon\test\vecmath　を試験ベンチとすべし
  改造をしたら、すくなくともその行の近傍に改造をした日時記述を入れておかないと
  外への悪影響があるときに発見しにくいよ。

  (C)Mathematical Assist Design Lab.
  (C)数理設計研究所　Hal.T
  1986.1    基本設計　関数タイプ
  1996/1/1  ANSI対応  クラステンプレート型
  1998/11/20 Hal.T  全体的な整備とANSI2への対応
  1999/10/5  Hal.T ラジアン、度のクラス整備

  2000/01
    Spehrical　は１引数では構築されず、３つ必要
    VecNsの参照宣言のときにめぐりめぐってdoubleではなくSphericalとして
    扱われ(0,0,0,1)になってしまったから　

  2004/12/15 WRITECSV　V2Dのバグ修正
  　
  //======================================================
  クラス名と内容
  //======================================================
  template  <class T> class Vec2;  //2要素ベクトル
  template  <class T> class Vec3;   //3要素ベクトル
  template  <class T> class VecN;   //N要素ベクトル ポインタ配列
  template  <class T, int jigen> class VecNs;  //N要素　静的配列
  template  <class T> class VecMtx; //N、M要素行列
  class   Spherical;                //球座標仮定義
  template  <class T>  class   SphericalVector //球座標＋別のデータ構造

  template  <class T>  int GetSign(T d) //符号判定 0を含む＋は整数１負側は-1
  int round(double x)  //	浮動小数点数の四捨五入による整数化
  inline double round(double x, int keta) //小数点以下の桁数指定
  template  <class T> double  Deg2Rad(T d) //度とラジアンの変換
  template  <class T> double  Rad2Deg(T d)
  template  <class T> class  DEG;         //角度扱い
  template  <class T> class  RADIAN
  //ファイル入出力
  template  <class T, template<class> class U> int ReadBin
  template  <class T, template<class> class U> int WriteBin
  //CSV入出力
  template  <class T, template<class> class U> int ReadCSV
  template  <class T, template<class> class U> int WriteCSV
  //======================================================

  ★２行２列以外では最初にIdentityで単位行列にする 2000/02/21 Hal.T
  void  RotateZ(double p)  など

  2,3,N次元ベクトル、球面座標の演算ライブラリ

  クラスの中に静的にベクトルを実現する場合にはコンストラクタに
  引数が持ち込めない。
  VecMtxなどでは次数初期化関数を用意してあるので、それで初期化する。

  コンパイラによりHファイル内にコードを全部書かなくてはいけないことがある
　そこでテンプレートらしくすべてを定義文としてＨファイルに納めた

　言葉の定義：
  学術図書出版　代数学と幾何学　P5
  「行列の横に並んでいる数の組を行といい、上から順に第一行、、、第ｍ行という。
  また横に並んでいる数の組を列といい、左から順に第一列、、第ｎ列という。
  したがってAijは第i行、第j列の成分または元素といい、行列自身は(Aij)のように
  略記する」

  Rowは行、columは列である（MathCAD,HELPより)　違う！
  行は column、　列はrow　である 2006/06/14 Hal.T

  関数テンプレートはクラス定義の中に書くと、クラスの内部関数としてしか
  使えないので、外に書く   abs など

--------使い方の注意-----------
★ポインタ型で配列を生成するとき
次数を持ち込めないので標準のコンストラクタをnewで生成した後に
  VecN.SetDimN(n)
  VecMtx.SetDimN(n,m)
で次元を確定させメモリを確保させる。

例： Ｎ次元ベクトル ４次のベクトルをｍ個ポインタで配列したいときには

  VecN<double>  *v31;    他の部分で宣言

  {
    if(v31)  delete [] v31;
    v31 = new VecN<double> [m];
    for(int n=0;  n<m;  v31[n++].SetDimN(4));
  }
★静的に配列を生成するとき
　Vec2<double> V2[256];
　ではヘッダレベルでなにも関数が実体化されていないのでエラーが起きる。
　ダミーとして
　Vec2<double> VP2, V2[256];
　とすれば生成される

★--------ライブラリ・メンテナンス上の注意-----------
class VecN<T> {
  VecN<T> operator+(VecN<T> &q)
  {  VecN<T> vt(*this);  for(int k=0;  k<vt.nt;  k++)  vt.vp[k]  += q.vp[k];
    return vt;  }
} a;
　これは
  d=a+b　では
  a.operator+(&b) として呼び出される

★ Vec2,Vec3,VecN,VecNs　
  operator=(T &src）がていぎされひとつの値で全要素が初期化される


========================================================
マニュアル部
角度の型： RADIAN と DEGがある
　RADIAN = DEG, DEG = RADIAN 相互の自動型変換と自動係数変換
  RADIAN *= K 定数倍
  RADIAN = RADIAN, RADIAN = int,double etc 代入
  adj() ０から２π未満、または０から３６０未満への調整

========================================================
------------------------------*/
#include  <values.h>
#include  <complex.h>
#include  <math.h>
//--------------------------------------
//StrictEqualが定義されていると"="による自動型変換が無効になりAssignのみになる
#define StrictEqual
//コンパイラオプション　VecNsに定数比較があるので、その警告を出させない
#pragma option -w-ccc
//--------------------------------------
//このファイルで定義されているもの
//定数
#define PI1     3.1415926535897932384626433832795  /* π   */

#ifndef PI
  #define PI  PI1
#endif

#define PI2     6.2831853071795864769252867665590  /* 2π  */
#define PI4     12.566370614359172953850573533118  /* 4π  */

#define PID2    1.5707963267948966192313216916397  /* π/2 */
#define PID18_  0.17453292519943295769236907684886 /* π/18 */
#define PID180  0.017453292519943295769236907684886 /* π/180 */
#define DEG2RAD PID180
#define RAD2DEG 57.295779513082320876798154814105  /* 180/π */
#define PIPI    9.8696044010893586188344909998761  /* π^2 */
//#pragma option  -Jgx
template  <class T> class Vec2;   //2要素ベクトル
// power,abs,~(=絶対値),+,-,*,/,==,!=,+=,-=,*=,/=,+p,-p
//[]による参照もできる

template  <class T> class Vec3;   //3要素ベクトル
//  Vec2への参照で初期化できる
// power,abs,~(=絶対値),+,-,*(スカラー積),^(ベクトル積),
// /,==,!=,+=,-=,*=,/=,+p,-p
//[]による参照もできる

//N要素ベクトル　静的配列
template  <class T, int jigen> class VecNs;
template  <class T> class VecN;   //N要素ベクトル ポインタ配列
// power,abs2,abs,~(=絶対値),+,-,*(スカラー積),^(ベクトル積),
// /,==,!=,+=,-=,*=,/=,+p,-p
// SetDimN()次元数の設定
// =(球座標からの変換), =Vec3よりの変換

template  <class T> class VecMtx; //N、M要素行列
// power,abs2,abs,~(=絶対値),+,-,*(スカラー積),^(ベクトル積),
// /,==,!=,+=,-=,*=,/=,+p,-p
// SetDimN()次元数の設定
//  ColumN, RowN  次数の取得
// Det,Det2,Det3　式の値
// Rotate 回転行列の作成 RotateV,X,Y,Z
// SolveC1,mul
// identity、単位行列の生成

class   Spherical;                //球座標仮定義
//--------------------------------------
//符号判定 0を含む＋側は整数で１負側は-1
template  <class T>  int GetSign(T d)  { return  (d>=0)? 1 : -1; }

//	浮動小数点数の四捨五入による整数化
int round(double x) { return x < 0.0 ? ceil( x - 0.5 ): floor( x + 0.5 ); }

//	浮動小数点数の四捨五入による整数化　小数点以下の桁数指定
inline double round(double x, int keta) {
  double  k = pow10(keta);
  x*=k;
  x = x < 0.0 ? ceil( x - 0.5 ): floor( x + 0.5 );
  return  x/=k;
}

//--------------------------------------
//度とラジアンの変換
template  <class T> double  Deg2Rad(T d)  { return  d*PID180;  }
template  <class T> double  Rad2Deg(T d)  { return  d*RAD2DEG;  }
//--------------------------------------
//角度扱い（度とラジアン）
#ifndef StrictEqual
template  <class T> class  DEG; //仮定義
template  <class T> class  RADIAN
{
  public:
  //data
  T a;

  // constructors
  RADIAN(T p) { a = p;  }
  RADIAN() {  a  = 0;  }

  //０から２πの範囲に収める
  void  adj(void)  {
    if(0<=a && a<PI2) return;
    int n = a/PI2;    if(a<0.0) n--;
    a -= n*PI2;
  }

  //四則演算
  RADIAN<T> operator+(RADIAN<T> &q) { return RADIAN<T>(a+q.a); }
  RADIAN<T> operator-(RADIAN<T> &q)	{ return RADIAN<T>(a-q.x); }
  RADIAN<T> operator*(T k)  { return RADIAN<T>(a*k);}  //スカラ倍
  RADIAN<T> operator/(T k)  { return RADIAN<T>(a/k);}  //スカラ除算

  int operator==(RADIAN<T> &p) { return p.a==a; }
  int operator!=(RADIAN<T> &p) { return p.a!=a; }

  //これ自身に対するオペレート
  RADIAN<T> &operator+=(RADIAN<T> &z) { a+=z.a; return *this; }
  RADIAN<T> &operator+=(T z) { a+=z; return *this; }
  RADIAN<T> &operator-=(RADIAN<T> &z) { a-=z.a; return *this; }
  RADIAN<T> &operator-=(T z) { a-=z; return *this; }
  RADIAN<T> &operator*=(T v)  { a *= v;  return *this; }
  RADIAN<T> &operator/=(T v)  { a /= v;  return *this; }

  RADIAN<T> operator+() { return *this; }
  RADIAN<T> operator-() { return RADIAN<T>(-x);  }

  RADIAN<T> Assign(DEG<T> &x) { a = Deg2Rad(x.a); return  *this;  }

    operator T() { return a; }  //自動型変換
    RADIAN<T> operator=(DEG<T> &x) {  return  Assign(x);  }
};

template  <class T> class  DEG
{
  public:
  //data
  T a;

  // constructors
  DEG(T p) { a = p;  }
  DEG() {  a  = 0;  }

  //０から360未満の範囲に収める
  void  adj(void)  {
    if(0<=a && a<360.0) return;
    int n = a/360.0;    if(a<0.0) n--;
    a -= n*360.0;
  }

  //四則演算
  DEG<T> operator+(DEG<T> &q) { return DEG<T>(a+q.a); }
  DEG<T> operator-(DEG<T> &q)	{ return DEG<T>(a-q.x); }
  DEG<T> operator*(T k)  { return DEG<T>(a*k);}  //スカラ倍
  DEG<T> operator/(T k)  { return DEG<T>(a/k);}  //スカラ除算

  int operator==(DEG<T> &p) { return p.a==a; }
  int operator==(DEG<T> p) { return p.a==a; }
  int operator!=(DEG<T> &p) { return p.a!=a; }
  int operator<(DEG<T> p) { return p.a<a; }

  //これ自身に対するオペレート
  DEG<T> &operator+=(T z) { a+=z; return *this; }
  DEG<T> &operator-=(T z) { a-=z; return *this; }
  DEG<T> &operator+=(DEG<T> &z) { a+=z.a; return *this; }
  DEG<T> &operator-=(DEG<T> &z) { a-=z.a; return *this; }
  DEG<T> &operator*=(T v)  { a *= v;  return *this; }
  DEG<T> &operator/=(T v)  { a /= v;  return *this; }

  DEG<T> operator+() { return *this; }
  DEG<T> operator-() { return DEG<T>(-x);  }

  DEG<T>  Assign(RADIAN<T> &x) { a = Rad2Deg(x.a); return  *this;  }

    operator T() { return a; }  //自動型変換
    DEG<T>  operator=(RADIAN<T> &x) { return  Assign(x);  }
};
#endif //StrictEqual   //}
#include  <fstream.h>
//--------------------------------------
//ファイル入出力
template  <class T, template<class> class U>
 int ReadBin(char *na,  U<T> *v, int n) {  return  v->ReadBin(na, v, n); }
template  <class T, template<class> class U>
 int WriteBin(char *na, U<T> *v, int n) {  return  v->WriteBin(na, v, n); }
//--------------------------------------
//CSV入出力
/*
  例：
  void __fastcall TForm1::Button25Click(TObject *Sender)
  {
  　V3D v3[3], vv[40];
　  v3[0].x=1;　v3[1].x=3;　v3[2].x=5;
  　WriteCSV("test.csv", v3, 3);
  　int n=ReadCSV("test.csv", vv, 4);

  VecNのCSV読み書きは動かない　VecNs　を使うこと　2005/05/27

  }
*/
// ファイル名、　読み取り場所、　読み取る数
template  <class T, template<class> class U>
  int ReadCSV(char *na,  U<T> *v, int n) {
  int m;
  ifstream  ifs(na, ios::in);  for(m=0; m<n &&  !ifs.eof();  m++)  ifs >> *v++;
  if(ifs.eof()) m--;
  ifs.close();
  return  m;  //読み取ったデータ数を返す
}

template  <class T, template<class> class U>
  int WriteCSV(char *na, U<T> *v, int n) {
  ofstream  ofs(na, ios::out);  for(; n>0;  n--)  ofs << *v++;
  ofs.close();
  return  n;  //書き残ったデータ数を返す
}
//--------------------------------------
/*固定小数点　16bit中央に小数点位置があるFFT用 intは32ビット
  デバッグ中だよお 2005/06/14
*/
class  FIX16_16 {
public:
  //data
  int x;
  // constructors
  FIX16_16() { x = 0;  }
  FIX16_16(double d) { x = d*65536; }
  FIX16_16(int  d) { x = d<<16; }
  FIX16_16(int  d,  unsigned int e) { x = d<<16 + e; }

  //絶対値
  FIX16_16  abs(void)   { return  x & 0x7FFFFFF; }

  //このクラスへの代入　自動変換
  FIX16_16  &operator=(const  int &d) {  x = d<<16;  return *this; }
  FIX16_16  &operator=(const  double  &d) { x = d*65536;  return  *this; }
  FIX16_16  &Assign(int d)  { x=d<<16;  return *this; }

  //このクラスからデフォルト型への代入　自動変換
  operator  int() const { return  x>>16; }
  operator  float() const { return  x/65536.0; }       //OK

  //四則演算  ここから下はまだ試験していない
  FIX16_16  operator+(FIX16_16  &q) { return FIX16_16(x+q.x); }

  FIX16_16  operator-(FIX16_16  &q)	{ return FIX16_16(x-q.x); }

  FIX16_16  operator*(FIX16_16  &q) { return  (x*q.x); }  //スカラー積

  FIX16_16  operator*(int k)  { return FIX16_16(x*k);}  //スカラ倍

  FIX16_16  operator/(int k)  { return FIX16_16(x/k);}  //スカラ除算

  //これ自身に対するオペレート
  FIX16_16  &operator+=(FIX16_16  &z) { x+=z.x; return *this; }
  FIX16_16  &operator-=(FIX16_16  &z) { x-=z.x; return *this; }
  FIX16_16  &operator*=(int v)  { x *= v;  return *this; }
  FIX16_16  &operator/=(int v)  { x /= v;  return *this; }

  FIX16_16  operator+() { return *this; }
  FIX16_16  operator-() { return FIX16_16(-x);  }
};
//--------------------------------------
//２次のベクトル
template  <class T>
class  Vec2 {
public:
  //data
  T x, y;

  // constructors
  Vec2(T p, T q) { x = p;  y = q;  }
  Vec2() {  x=y=0;  }

  //パワー値
  T power(void) { return  (x*x + y*y); }
  T power(Vec2<T>  &v) { return  (v.x*v.x + v.y*v.y); }

  //絶対値
  T operator~() { return  sqrt(x*x + y*y); }
  T abs(void)   { return  sqrt(x*x + y*y); }
  T abs(Vec2<T>  &v) { return  sqrt(v.x*v.x + v.y*v.y); }

  //規格化　長さを１にする
  void  Standard(void)  { *this /= abs(*this);  }


  //四則演算
  Vec2<T> operator+(Vec2<T> &q) { return Vec2<T>(x+q.x, y+q.y); }

  Vec2<T> operator-(Vec2<T> &q)	{ return Vec2<T>(x-q.x, y-q.y); }

  T operator*(Vec2<T> &q) { return  (x*q.x + y*q.y); }  //スカラー積

  Vec2<T> operator*(T k)  { return Vec2<T>(x*k, y*k);}  //スカラ倍

  Vec2<T> operator/(T k)  { return Vec2<T>(x/k, y/k);}  //スカラ除算

  //int operator==(Vec2<T> &p) { return p.x==x && p.y==y; }
  //int operator!=(Vec2<T> &p)  { return p.x!=x || p.y!=y; }

  //ベクトル積で面積、3DのZ値を帰す
  double operator^(Vec2<T> &q)  { return  x*q.y-y*q.x; }

  //これ自身に対するオペレート
  Vec2<T> &operator+=(Vec2<T> &z) { x+=z.x; y+=z.y; return *this; }
  Vec2<T> &operator-=(Vec2<T> &z) { x-=z.x; y-=z.y; return *this; }
  Vec2<T> &operator*=(T v)  { x *= v;  y*=v;  return *this; }
  Vec2<T> &operator/=(T v)  { x /= v;  y/=v;  return *this; }

  Vec2<T> operator+() { return *this; }
  Vec2<T> operator-() { return Vec2<T>(-x, -y);  }

  //x y の代わりの[]引数による参照
  T& operator[](int n)  { return  (n==0)? x : y;  }

  //代入　
  Vec2<T> &Assign(T &v)  { x=y=v;  return *this; } //全て同じ値
  Vec2<T> &Assign(Vec2<T> &v)  { x=v.x;  y=v.y;  return *this; }
  Vec2<T> &Assign(Vec3<T> &v)  { x=v.x;  y=v.y;  return *this; }
  Vec2<T> &Assign(VecNs<T,4> &v) { x=v[0];  y=v[1];  return *this; }
  Vec2<T> &Assign(VecN<T> &v)  { x=v[0]; y=v[1]; return *this; }

  friend  ostream &operator<<(ostream &s, Vec2<T> &p) {
    s << p.x << ", " << p.y << endl; return  s;
  }

  friend  istream &operator>>(istream &s, Vec2<T> &p) {
    char buf[1024];         double  x, y;
    s.getline(buf, 1024);   sscanf(buf, "%lf,%lf", &x , &y);
    p.x = x;    p.y = y;
    return  s;
  }

  int WriteBin(FILE *fp) {  return  fwrite(&x,  sizeof(T), 2, fp);  }
  int WriteBin(FILE *fp,  Vec2<T> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  Vec2<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {  return  fread(&x,  sizeof(T), 2, fp);  }
  int ReadBin(FILE *fp,  Vec2<T> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  Vec2<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

};
/*
  関数テンプレートはクラス定義の中に書くと、クラスの内部関数としてしか
  使えないので、外に書く
*/

template <class T> T abs(Vec2<T>  v) { return  sqrt(v.x*v.x + v.y*v.y); };
template <class T> T power(Vec2<T>  v) { return  (v.x*v.x + v.y*v.y); }

//STLのuniqueを使うために関数テンプレートに変えた 2001/2/21 Hal.T 
template  <class T> bool operator==(const Vec2<T> &p, const Vec2<T> &q)  {
  return p.x==q.x && p.y==q.y; }
template  <class T> bool operator!=(const Vec2<T> &p, const Vec2<T> &q)  {
  return p.x!=q.x || p.y!=q.y; }
//--------------------------------------
//3次のベクトル
//ベクトル積などの計算順序は使うときに（）で保証すること
template  <class T>
class  Vec3 {
  public:
  //data
  T x, y, z;

  // constructors
  Vec3() {  x=y=z=0;  }
  Vec3(T p, T q, T r) { x=p;  y=q;  z=r;  }
  Vec3(Vec2<T> &v2) { x=v2.x;  y=v2.y;  z=0.0;  }

  //パワー値
  T power(void) { return  (x*x + y*y + z*z); }
  T power(Vec3<T>  &v) { return  (v.x*v.x + v.y*v.y + v.z*v.z); }

  //絶対値
  T operator~() { return  sqrt(x*x + y*y + z*z); }
  T abs(void) const  { return  sqrt(x*x + y*y + z*z); }
  T abs(const Vec3<T>  &v)  {  return  sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }

  //規格化　長さを１にする
  void  Standard(void)  { *this /= abs(*this);  }

  //四則演算
  Vec3<T> operator+(Vec3<T> &q) { return Vec3<T>(x+q.x, y+q.y,  z+q.z); }

  Vec3<T> operator-(Vec3<T> &q)	{ return Vec3<T>(x-q.x, y-q.y,  z-q.z); }

  T operator*(Vec3<T>  &q)  { return  (x*q.x + y*q.y + z*q.z); }  //スカラー積

  Vec3<T> operator^(Vec3<T> &q) //ベクトル積
    { return  Vec3<T>(y*q.z-z*q.y, z*q.x-x*q.z, x*q.y-y*q.x); }

  Vec3<T> operator*(T k)  { return Vec3<T>(x*k, y*k,  z*k); }     //スカラ倍

  Vec3<T> operator/(T k)  { return Vec3<T>(x/k, y/k, z/k); }   //スカラー除算

  //bool operator==(const Vec3<T> &p)  { return p.x==x && p.y==y && p.z==z; }
  //bool operator!=(const Vec3<T> &p)  { return p.x!=x || p.y!=y || p.z!=z; }

  //これ自身に対するオペレート
  Vec3<T> &operator+=(Vec3<T> &a) { x+=a.x; y+=a.y; z+=a.z; return *this; }
  Vec3<T> &operator-=(Vec3<T> &a) { x-=a.x; y-=a.y; z-=a.z; return *this; }
  Vec3<T> &operator*=(T v)  { x*=v;  y*=v;  z*=v; return *this; }
  Vec3<T> &operator/=(T v)  { x/=v;  y/=v;  z/=v; return *this; }

  Vec3<T> operator+() { return *this; }
  Vec3<T> operator-() { return Vec3<T>(-x, -y,  -z);  }

  Vec3<T>& operator=(const  Vec3<T> &p) {
    x=p.x;  y=p.y;  z=p.z;  return *this; }

  Vec3<T>& round(int keta)
    { x=::round(x, keta); y=::round(y,keta); z=::round(z,keta); return *this; }

  Vec3<T>& Reciprocal(void)  { x=1/x; y=1/y, z=1/z; return *this; } //逆数化
  Vec3<T>& Reciprocal(const T &src)
    {  x=src/x; y=src/y;  z=src/z; return *this;  } //逆数化 src/x

  //型の違うものの代入はAssignとする
  Vec3<T>&  Assign(const T &src) { x=src;  y=src;  z=src;    return *this; }
  Vec3<T>&  Assign(const Vec2<T> &p) { x=p.x;  y=p.y;  z=0;  return *this; }
  Vec3<T>&  Assign(const Vec3<T> &p) { x=p.x;  y=p.y;  z=p.z;  return *this; }
  Vec3<T>&  Assign(const VecN<double> &p) { x=p.vp[0];  y=p.vp[1];  z=p.vp[2];  return  *this;  }
  Vec3<T>&  Assign(const  VecNs<double,3> &src) {
    x=src.v[0];  y=src.v[1];  z=src.v[2];    return *this; }
  Vec3<T>&  Assign(const  VecNs<double,4> &src) {
    x=src.v[0];  y=src.v[1];  z=src.v[2];    return *this; }
  //球座標から直交座標に変換する
  Vec3<T>&  Assign(const Spherical &src)  {
    double  rxy;
    x=(rxy=src.r*cos(src.theta))*cos(src.phi);    y=rxy*sin(src.phi);
    z=src.r*sin(src.theta);   return *this;
  }

  #ifndef StrictEqual  //{
  Vec3<T>& operator=(const  T &src) { return  Assign(src);  }
//  Vec3<T>& operator=(const  Vec3<double> &src) {  return  Assign(src);  }
//  Vec3<T>& operator=(const  Vec3<int> &src) { return  Assign(src);  }
  Vec3<T>& operator=(const  VecN<double> &src) { return  Assign(src);  }
  Vec3<T>& operator=(const  VecNs<double,3> &src) { return  Assign(src);  }
  Vec3<T>& operator=(const  VecNs<double,4> &src) { return  Assign(src);  }
  Vec3<T>& operator=(const Spherical &src)  { return  Assign(src);  }
  #endif //}

  T& operator[](int n)  { return  (n==0)? x : ((n==1)? y : z);  }

  friend  ostream &operator<<(ostream &s, Vec3<T> &p) {
    s << p.x << ", " << p.y << ", " << p.z << endl; return  s;
  }

  friend  istream &operator>>(istream &s, Vec3<T> &p) {
    char buf[1024];         double  x, y, z;
    s.getline(buf, 1024);   sscanf(buf, "%lf,%lf,%lf", &x , &y, &z);
    p.x = x;    p.y = y;    p.z = z;
    return  s;
  }

  int WriteBin(FILE *fp) {  return  fwrite(&x,  sizeof(T), 3, fp);  }
  int WriteBin(FILE *fp,  Vec3<T> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  Vec3<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {  return  fread(&x,  sizeof(T), 3, fp);  }
  int ReadBin(FILE *fp,  Vec3<T> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  Vec3<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }
};

template  <class T> T power(Vec3<T>  v) { return  (v.x*v.x + v.y*v.y + v.z*v.z); }
template  <class T> T abs(Vec3<T>  v) {  return  sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }

//STLのuniqueを使うために関数テンプレートに変えた 2001/2/21 Hal.T 
template  <class T> bool operator==(const Vec3<T> &p, const Vec3<T> &q)  {
  return p.x==q.x && p.y==q.y && p.z==q.z; }
template  <class T> bool operator!=(const Vec3<T> &p, const Vec3<T> &q)  {
  return p.x!=q.x || p.y!=q.y || p.z!=q.z; }


//--------------------------------------
/*Ｎ次元ベクトル ポインタ配列型
  ４次のベクトルをｍ個ポインタで配列したいときには

  VecN<double>  *v31;    他の部分で宣言

  {
    if(v31)  delete [] v31;
    v31 = new VecN<double> [m];
    for(int n=0;  n<m;  v31[n++].SetDimN(4));
  }
*/
template  <class T> class VecN {
public:
  int vn; //要素数
  T  *vp; //ベクトルデータへのポインタ

  // constructors
  VecN()  { vp=0; vn=0; }

  VecN(const VecN<T> &v)  //既存の内容で初期化する
  {  vp = new T[vn = v.vn];  for(int k=0;  k<vn; k++) vp[k] =  v.vp[k];}

  VecN(int n, T *p) //既存の内容で初期化する
  {  vp = new T[vn = n];  if(p) { for(int k=0;  k<n; k++) vp[k] =  p[k];  }}

  VecN(int n, T d=0)  //引数で初期化する
  {  vp = new T[vn = n];  for(int k=0;  k<vn; k++) vp[k] =  d;}

  ~VecN() { DeleteVecN(); }

  //関数
  void  DeleteVecN(void) { if(vp) delete [] vp; vp=0;};
  //次元数を設定
  void  SetDim(int  n) { SetDimN(n);  }
  void  SetDimN(int  n) { DeleteVecN();  vp = new T  [vn = n];  }

  int Dim(void) { return  vn; };                     //次元（要素数）を返す

  friend  int GetDim(VecN<T> &p) { return p.vn;  }   //次元（要素数）を返す

  void  Clear(T d)  { //全部の初期化
    for(int k=0;  k<vn; k++) vp[k] =  d;
  }

  //絶対値
  T operator~()
  { T d=0.0; for(int k=0; k<vn; k++) d+=(T)(vp[k]*vp[k]); return sqrt(d);}
  T abs(void)
  { T d=0.0; for(int k=0; k<vn; k++) d+=(T)(vp[k]*vp[k]); return sqrt(d);}
  T abs(VecN<T> &v)
  { T d=0.0; for(int k=0; k<v.vn; k++) d+=(T)(v.vp[k]*v.vp[k]); return sqrt(d);}

  //絶対値で開平していないもの
  T power(void)
  { T d=0.0;  for(int k=0;  k<vn; k++)  d+=(T)(vp[k]*vp[k]); return  d; }
  T power(VecN<T> &v)
  { T d=0.0;  for(int k=0;  k<v.vn; k++)  d+=(T)(v.vp[k]*v.vp[k]); return  d; }

  //規格化　長さを１にする
  void  Standard(void)  { *this /= abs(*this);  }

  VecN<T> operator+(VecN<T> &q)
  { VecN<T> vt(*this); for(int k=0; k<vt.vn; k++) vt.vp[k]+=q.vp[k]; return vt;}

  VecN<T> operator-(VecN<T> &q)
  { VecN<T> vt(*this); for(int k=0; k<vt.vn; k++) vt.vp[k]-=q.vp[k]; return vt;}

  T operator*(VecN<T>  &q)  //スカラ積
  { T t=0; for(int k=0; k<vn; k++) t+=vp[k]*q.vp[k]; return t;}

  VecN<T> operator*(T a)
  { VecN<T> vt(*this); for(int k=0; k<vn; k++) vt.vp[k]*=a; return vt;}

  VecN<T> operator/(T a)
  { VecN<T> vt(*this); for(int k=0; k<vn; k++) vt.vp[k]/=a; return vt;}

  T& operator[](int n)  { return  (*this).vp[n];  }

  VecN<T> &operator+=(VecN<T> const &p)
  { for(int k=0;  k<vn;  k++)  vp[k]  += p.vp[k]; return *this;  }

  VecN<T> &operator-=(VecN<T> const &p)
  { for(int k=0;  k<vn;  k++)  vp[k]-=p.vp[k]; return *this;  }

  VecN<T> &operator*=(T a)
  { for(int k=0;  k<vn;  k++)  vp[k]*=a;    return *this;  }

  VecN<T> &operator/=(T a)
  { for(int k=0;  k<vn;  k++)  vp[k]/=a;    return *this;  }


  VecN<T> operator+() { return  *this;  }   //＋符号
  VecN<T> operator-()                    //−符号
  {  VecN<T> vt(this->vn);  for(int k=0; k<vn; k++) vt[k]=-vp[k]; return vt; }

  VecN<T>& operator=(const T &src)
  { for(int k=0;  k<vn;  k++)  vp[k]=src; return *this;  }

  VecN<T>& operator=(const VecN<T> &src)
  { for(int k=0;  k<vn;  k++)  vp[k]=src.vp[k]; return *this;  }

  VecN<T>& operator=(const Vec3<T> &src)
  { vp[0]=src.x;  vp[1]=src.y;  vp[2]=src.z;    return *this;  }


  //納まる範囲でコピーする
  VecN<T>&  Assign(const  Vec2<T> &v2)  {
    vp[0]=v2.x;  if(vn>=2) vp[1] = v2.y;  return  *this;  }
  VecN<T>&  Assign(const  Vec3<T> &v3)  {
    vp[0]=v3.x;  if(vn>=2) vp[1] = v3.y;  if(vn>=3) vp[2] = v3.z;
    return  *this;
  }

  VecN<T>&  Assign(const  VecN<T> &vv)  {
    int n; for(n=0;  n<vv.vn && n<vn;  n++)  vp[n] = vv[n];
    return  *this;
  }

  //球面座標からの変換　アフィン変換用  最後の項は１にする
  #ifndef StrictEqual  //{
  VecN<T>& operator=(const Spherical &src)  { return  Assign(src);  }
  #endif  //}

  VecN<T>&  Assign(const Spherical &src)  {
    double  rxy;
    vp[0]=(rxy=src.r*cos(src.theta))*cos(src.phi);
    vp[1]=rxy*sin(src.phi);
    vp[2]=src.r*sin(src.theta);
    if(vn>3) vp[3]=1.0;
    return *this;
  }

  friend  ostream &operator<<(ostream &s, VecN<T> &p){  //2005/05/27　BUG不可解
    s << p.vp[0];
    for(int n=1;  n<p.vn; ) {  s << ", " << p.vp[n++];  }
    s << endl;  return  s;
  }

  friend  istream &operator>>(istream &s, VecN<T> &p) {
    char *cp, buf[1024];         double  x;
    s.getline(buf, 1024);
    for(int n=0, cp=buf;  n<p.vn; ) {
      p.v[n++]=atof(cp);   cp=strchr(cp, ',');
      if(cp)  cp++;
    }
    return  s;
  }

  //以下の二つの関数は基本レベルの変数タイプだけを扱う char int double

  int WriteBin(FILE *fp) {  return  fwrite(vp,  sizeof(T), vn, fp);  }
  int WriteBin(FILE *fp,  VecN<T> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  VecN<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {  return  fread(vp,  sizeof(T), vn, fp);  }
  int ReadBin(FILE *fp,  VecN<T> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  VecN<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

};

template  <class T>  T abs(VecN<T> v)
  { T d=0.0; for(int k=0; k<v.vn; k++) d+=(T)(v.vp[k]*v.vp[k]); return sqrt(d);}
template  <class T>  T power(VecN<T> v)
  { T d=0.0;  for(int k=0;  k<v.vn; k++)  d+=(T)(v.vp[k]*v.vp[k]); return  d; }

//STLのuniqueを使うために関数テンプレートに変えた 2001/2/21 Hal.T
template  <class T> bool operator==(const VecN<T> &p, const VecN<T> &q)  {
  int k=0, t=-1; for(; k<p.vn && t; k++) { t=(p.vp[k]==q.vp[k]); } return t;}
template  <class T> bool operator!=(const VecN<T> &p, const VecN<T> &q)  {
   for(int k=0; k<p.vn; k++) { if(p.vp[k]!=q.vp[k]) return -1; } return 0;}
//--------------------------------------
/*Ｎ次元ベクトル コンパイル時メモリ確保の型

構築方法
  VecNs<double,4> v4()
  VecNs<double,4> v4(v4a)     他への参照
  VecNs<double,4> v4=v4a      #
  VecNs<double,4> v4(3.3)     数値データ

  double  p[4] = {1,2,3,4};
  VecNs<double,4> v4(p)       数値配列へのポインタ

*/
template  <class T, int jigen> class VecNs {
  public:
  T v[jigen]; //ベクトルデータへのポインタ

  // constructors
  VecNs()  {  for(int n=0;  n<jigen; v[n++]=0); }
  //初期化データを持ちこむ
  VecNs(T d)  { for(int n=0;  n<jigen; v[n++]=d); }
  VecNs(T *d)  { for(int n=0;  n<jigen; v[n++]=*d++); }

  //関数
  int Dim(void) { return  jigen; };                     //次元（要素数）を返す

  //絶対値
  T operator~()
  { T d=0.0; for(int k=0; k<jigen; k++) d+=(T)(v[k]*v[k]); return sqrt(d);}
  T abs(void)
  { T d=0.0; for(int k=0; k<jigen; k++) d+=(T)(v[k]*v[k]); return sqrt(d);}
  T abs(VecNs<T, jigen> &v)
  { T d=0.0; for(int k=0; k<jigen; k++) d+=(T)(v.v[k]*v.v[k]); return sqrt(d);}

  //絶対値で開平していないもの
  T power(void)
  { T d=0.0;  for(int k=0;  k<jigen; k++)  d+=(T)(v[k]*v[k]); return  d; }
  T power(VecNs<T, jigen> &v)
  { T d=0.0;  for(int k=0;  k<jigen; k++)  d+=(T)(v.v[k]*v.v[k]); return  d; }

  //規格化　長さを１にする
  void  Standard(void)  { *this /= abs(*this);  }

  VecNs<T, jigen> operator+(VecNs<T, jigen> &q)
  { VecNs<T, jigen> vt(*this); for(int k=0; k<jigen; k++) vt.v[k]+=q.v[k]; return vt;}

  VecNs<T, jigen> operator-(VecNs<T, jigen> &q)
  { VecNs<T, jigen> vt(*this); for(int k=0; k<jigen; k++) vt.v[k]-=q.v[k]; return vt;}

  T operator*(VecNs<T, jigen>  &q)  //スカラ積
  { T t=0; for(int k=0; k<jigen; k++) t+=v[k]*q.v[k]; return t;}

  VecNs<T, jigen> operator*(T a)
  { VecNs<T, jigen> vt(*this); for(int k=0; k<jigen; k++) vt.v[k]*=a; return vt;}

  VecNs<T, jigen> operator/(T a)
  { VecNs<T, jigen> vt(*this); for(int k=0; k<jigen; k++) vt.v[k]/=a; return vt;}
/*
  int operator==(VecNs<T, jigen> &q)
  { int k=0, t=-1; for(; k<jigen && t; k++) { t=(v[k]==q.v[k]); } return t;}

  int operator!=(VecNs<T, jigen> &q)
  { for(int k=0; k<jigen; k++) { if(v[k]!=q.v[k]) return -1; } return 0;}
*/
  T& operator[](int n)  { return  v[n];  }

  VecNs<T, jigen> &operator+=(VecNs<T, jigen> const &p)
  { for(int k=0;  k<jigen;  k++)  v[k]  += p.v[k]; return *this;  }

  VecNs<T, jigen> &operator-=(VecNs<T, jigen> const &p)
  { for(int k=0;  k<jigen;  k++)  v[k]-=p.v[k]; return *this;  }

  VecNs<T, jigen> &operator*=(T a)
  { for(int k=0;  k<jigen;  k++)  v[k]*=a;    return *this;  }

  VecNs<T, jigen> &operator/=(T a)
  { for(int k=0;  k<jigen;  k++)  v[k]/=a;    return *this;  }


  VecNs<T, jigen> operator+() { return  *this;  }   //＋符号
  VecNs<T, jigen> operator-()                    //−符号
  {  VecNs<T, jigen> vt;  for(int k=0; k<jigen; k++) vt[k]=-v[k]; return vt; }

  VecNs<T, jigen>& operator=(const VecNs<T, jigen> &src)
  { for(int k=0;  k<jigen;  k++)  v[k]=src.v[k]; return *this;  }

  //同一型による要素初期化　行列作成に必要
  VecNs<T, jigen>& operator=(const T &src)
  { for(int k=0;  k<jigen;  k++)  v[k]=src; return *this;  }

  #ifndef StrictEqual //{
  VecNs<T, jigen>& operator=(const Vec3<T> &src)
  { v[0]=src.x;  v[1]=src.y;  v[2]=src.z;    return *this;  }

  //球面座標からの変換　アフィン変換用  最後の項は１にする
  VecNs<T,  jigen>& operator=(const Spherical &src)  { return  Assign(src);  }
  #endif  //}

  //球座標からのコピー　アフィン変換用に４次なら４次目は１にする　
  VecNs<T,  jigen>&  Assign(const Spherical &src)  {
    double  rxy;
    v[0]=(rxy=src.r*cos(src.theta))*cos(src.phi);    v[1]=rxy*sin(src.phi);
    v[2]=src.r*sin(src.theta);                       if(jigen>3) v[3]=1.0;
    return *this;
  }

  VecNs<T,  jigen>&  Assign(const  Vec2<T> &v2)  {
    v[0] = v2.x;  if(jigen>=2) v[1] = v2.y;  return *this;  }

  VecNs<T,  jigen>&  Assign(const  Vec3<T> &v3) {
    v[0] = v3.x;  if(jigen>=2) v[1] = v3.y;  if(jigen>=3) v[2] = v3.z;
    return  *this;
  }

  VecNs<T,  jigen>&  Assign(const  VecNs<T, jigen> &vv) { //同じサイズのコピー
    for(int n=0;  n<jigen;  n++)  v[n] = vv.v[n];
    return  *this;
  }

  //配列からのコピー 2000/02/13 Hal.T
  VecNs<T,  jigen>&  Assign(const  T  *p) {
    for(int n=0;  n<jigen;  n++)  v[n] = p[n];
    return  *this;
  }

  friend  ostream &operator<<(ostream &s, VecNs<T, jigen> &p) {
    s << p.v[0];
    for(int n=1;  n<jigen; ) {  s << ", " << p.v[n++];  }
    s << endl;  return  s;
  }

  friend  istream &operator>>(istream &s, VecNs<T, jigen> &p) {
    char *cp, buf[1024];         double  x;
    s.getline(buf, 1024);   cp = buf;
    for(int n=0; n<jigen && cp; ) {
      p.v[n++]=atof(cp);   cp=strchr(cp, ',');
      if(cp)  cp++;
    }
    return  s;
  }


  //以下の二つの関数は基本レベルの変数タイプだけを扱う char int double

  int WriteBin(FILE *fp) {  return  fwrite(v,  sizeof(T), jigen, fp);  }
  int WriteBin(FILE *fp,  VecNs<T,  jigen> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  VecNs<T,  jigen> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {  return  fread(v,  sizeof(T), jigen, fp);  }
  int ReadBin(FILE *fp,  VecNs<T,  jigen> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  VecNs<T,  jigen> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

};
template  <class T, int jigen>
  int ReadCSVs(char *na,  VecNs<T, jigen> *v, int n) {
  int m;
  ifstream  ifs(na, ios::in);  for(m=0; m<n &&  !ifs.eof();  m++)  ifs >> *v++;
  if(ifs.eof()) m--;
  ifs.close();
  return  m;  //読み取ったデータ数を返す
}
template  <class T, int jigen>
  int WriteCSVs(char *na,  VecNs<T, jigen> *v, int n) {
  ofstream  ofs(na, ios::out);  for(; n>0;  n--)  ofs << *v++;
  ofs.close();
  return  n;  //書き残ったデータ数を返す
}

//絶対値
template  <class T, int jigen>  T abs(VecNs<T, jigen> &v)
{ T d=0.0; for(int k=0; k<jigen; k++) d+=(T)(v.v[k]*v.v[k]); return sqrt(d);}

//絶対値の自乗
template  <class T, int jigen>  T power(VecNs<T, jigen> &v)
{ T d=0.0;  for(int k=0;  k<jigen; k++)  d+=(T)(v.v[k]*v.v[k]); return  d; }

//関数テンプレート　ポインタへのコピー  2000/02/13 Hal.T
template  <class T, int jigen>  T * AssignCopy(T *p, VecNs<T, jigen> &v)
{  T *q=p; for(int k=0;  k<jigen; k++)  *q++=(T)(v.v[k]);  return  p; }

//STLのuniqueを使うために関数テンプレートに変えた 2001/2/21 Hal.T
template  <class T, int jigen> bool operator==(const VecNs<T, jigen> &p, const VecNs<T, jigen> &q)  {
  int k=0, t=-1; for(; k<jigen && t; k++) { t=(p.vp[k]==q.vp[k]); } return t;}
template  <class T, int jigen> bool operator!=(const VecNs<T, jigen> &p, const VecNs<T, jigen> &q)  {
   for(int k=0; k<jigen; k++) { if(p.vp[k]!=q.vp[k]) return -1; } return 0;}

//--------------------------------------
//VecMtx n行 m列　の直方行列 VecN を使った行列
template  <class T>
class  VecMtx  //: public  VecN<T>
{
  private:
  public:
  int Row, Colum;
  VecN<T> *vp;

  // constructors
  VecMtx()  { Colum=Row=0;  vp=0; }  //次数の定義 SetDim(c,r)を呼び出すこと

  VecMtx(int row,  int  colum){
    Colum = colum;    Row = row;  vp  = new VecN<T> [Row];
    for(int k = 0;  k<Row;  k++)  vp[k].SetDimN(Colum); }

  VecMtx(const VecMtx &m){
    Colum = m.Colum;    Row = m.Row;  vp  = new VecN<T> [Row];
    for(int k = 0;  k<Row;  k++)  {  vp[k].SetDimN(Colum); vp[k] = m.vp[k];  }}

  ~VecMtx() {   delete  []  vp; }

  //次数の(再)定義
  void  DeleteVecN(void) { if(vp!=0) delete [] vp; vp=0;};

  void  SetDim(int  row,  int colum) { SetDimN(row, colum);  }
  void  SetDimN(int  row, int colum) {
    DeleteVecN();    Colum = colum;    Row = row;
    vp  = new VecN<T> [Row];
    for(int k = 0;  k<Row;  k++)  vp[k].SetDimN(Colum);
  }

  int RowN(void) { return Row;  }     //行数
  int ColumN(void) { return Colum;  } //列数

  void  Clear(T d)  { //全部の初期化
    for(int k = 0;  k<Row;  k++)  vp[k].Clear(d);
  }

  void  Identity(void)  { //単位行列の生成
    if(Colum != Row)  return;
    for(int n=0;  n<Row;  n++)  { for(int m=0;  m<Colum;  m++)  vp[n][m]=(T)0; }
    for(int n=0;  n<Row;  n++)  vp[n][n]=1;
  }

  //行列式の値、２か３の正方行列　それ以外の場合は最大値
  double Det2(void) {
    return  vp[0][0]*vp[1][1]-vp[0][1]*vp[1][0];
  }

  double Det3(void) {
    return  vp[0][0]*vp[1][1]*vp[2][2]
            +vp[0][1]*vp[1][2]*vp[2][0]
            +vp[0][2]*vp[1][0]*vp[2][1]
            -vp[0][2]*vp[1][1]*vp[2][0]
            -vp[0][1]*vp[1][0]*vp[2][2]
            -vp[0][0]*vp[1][2]*vp[2][1];

  }

  double Det4(void) {
    return
     vp[0][0]*(vp[1][1]*(vp[2][2]*vp[3][3]-vp[2][3]*vp[3][2])
              +vp[1][2]*(vp[2][3]*vp[3][1]-vp[2][1]*vp[3][3])
              +vp[1][3]*(vp[2][1]*vp[3][2]-vp[2][2]*vp[3][1]))
    -vp[0][1]*(vp[1][0]*(vp[2][2]*vp[3][3]-vp[2][3]*vp[3][2])
              +vp[1][2]*(vp[2][3]*vp[3][0]-vp[2][0]*vp[3][3])
              +vp[1][3]*(vp[2][0]*vp[3][2]-vp[2][2]*vp[3][0]))
    +vp[0][2]*(vp[1][0]*(vp[2][1]*vp[3][3]-vp[2][3]*vp[3][1])
              +vp[1][1]*(vp[2][3]*vp[3][0]-vp[2][0]*vp[3][3])
              +vp[1][3]*(vp[2][0]*vp[3][1]-vp[2][1]*vp[3][0]))
    -vp[0][3]*(vp[1][0]*(vp[2][1]*vp[3][2]-vp[2][2]*vp[3][1])
              +vp[1][1]*(vp[2][2]*vp[3][0]-vp[2][0]*vp[3][2])
              +vp[1][2]*(vp[2][0]*vp[3][1]-vp[2][1]*vp[3][0]));
  }

  double Det(void) {
    if(Colum==2 && Row==2)  return  Det2();
    else  if(Colum==3 && Row==3)  return  Det3();
    else  if(Colum==4 && Row==4)  return  Det4();
    else  return  MAXDOUBLE;
  }

  //２次元回転行列
  void  RotateV(double p)  {
    vp[0][0] = vp[1][1] = cos(p);
    vp[1][0] = -(vp[0][1] = sin(p));
  }

  //任意の回転ベクトルvを軸にする(p)回転行列要素の生成
  //日刊工業新聞　「入門　ＣＡＤ図形計算手法」ｐ３２
  //３行３列以外では最初にIdentityで単位行列にしてから実行する
  void  RotateV(Vec3<T> &v3, double p)  {
    double  c, c1, s;
    Vec3<T> v = v3;
    c1 = 1.0-(c = cos(p));  s = sin(p);
    if(v.power()==0.0)  {  Identity(); return; }
    v.Standard();
    vp[0][0]=c1*v.x*v.x+c;     vp[0][1]=c1*v.x*v.y-v.z*s; vp[0][2]=c1*v.x*v.z+v.y*s;
    vp[1][0]=c1*v.x*v.y+v.z*s; vp[1][1]=c1*v.y*v.y+c;     vp[1][2]=c1*v.y*v.z-v.x*s;
    vp[2][0]=c1*v.x*v.z-v.y*s; vp[2][1]=c1*v.y*v.z+v.x*s; vp[2][2]=c1*v.z*v.z+c;
  }

  void  RotateVA(Vec3<T> &v3)  {  RotateV(v3, abs(v3));  }

  //回転表現 X,Y,Z　要素だけを変化させる　
  //３行３列以外では最初にIdentityで単位行列にしてから実行する
  void  RotateZ(double  x,  double  y)  {
    vp[0][0] = vp[1][1] = x;    vp[1][0] = -(vp[0][1]= y);
  }
  //２行２列以外では最初にIdentityで単位行列にしてから実行する 2000/02/21 Hal.T
  void  RotateZ(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    vp[0][0] = vp[1][1] = cos(p);    vp[1][0] = -(vp[0][1]= sin(p));
  }
  void  RotateX(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    vp[1][1] = vp[2][2] = cos(p);    vp[1][2] = -(vp[2][1]= sin(p));
  }
  void  RotateY(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    vp[0][0] = vp[2][2] = cos(p);    vp[0][2] = -(vp[2][0]= sin(p));
  }

  void  SolveC1(void)  // 対角を１にして、係数を第一列に吐き出す
  {  int nC1, nC2;  T  tmpCC;
    for(nC1 = 1;  nC1 < Row; nC1++) {
      tmpCC = vp[nC1][nC1];
      for(nC2 = 1;  nC2 < Row;  nC2++)  {
        if(nC1 != nC2)  vp[nC2] -=  vp[nC1]*(vp[nC2][nC1]/tmpCC);
      }
    }

    for(nC1=1;  nC1 < Row;  nC1++)  { // 行列の対角要素を１に規格化する
      tmpCC = 1.0/vp[nC1][nC1];
      vp[nC1] *=  tmpCC;
    }
  }

  VecN<T>& operator[](int n)  { return  (*this).vp[n];  }

  //Tの内部仕様にもよるが単純コピーでよければColum*Rowの量だけmemcpyすればよい
  VecMtx<T>& operator=(const  VecMtx<T> &src)
  {
    if(src.Colum != Colum || src.Row != Row) return  *this;
    for(int k = 0;  k<Row;  k++)  vp[k] = src.vp[k];
    return *this;
  }

  //定数倍
//  friend VecMtx<T> operator*(VecMtx<T> &a, T d){  VecMtx<T> c(a.Row, a.Colum);
//  for(int m=0;  m<c.Row;  m++)  {    c[m] = a[m]*d;  }  return c;}
  VecMtx<T> operator*(T d)  {
    VecMtx<T> &a = *this;
    VecMtx<T> c(a.Row, a.Colum);
    for(int m=0;  m<c.Row;  m++)  {    c[m] = a[m]*d;  }  return c;
  }


  //VecMtx*VecMtx　行列の乗算
  //行の次数があわなければ零で埋まった行列を帰す、内部で関数版を使用している
  //代数学と幾何学（学術図書出版)P11
  VecMtx<T> operator*(VecMtx<T> &b)
  { VecMtx<T> &a = *this;
    VecMtx<T> c(a.Row, a.Colum);    c.mul(a,b);  return c;}

  Vec3<T> operator*(Vec3<T> &b)
  {
    VecMtx<T> &a = *this;
    Vec3<T> c;
    c.x = a[0][0]*b.x+a[0][1]*b.y+a[0][2]*b.z;
    c.y = a[1][0]*b.x+a[1][1]*b.y+a[1][2]*b.z;
    c.z = a[2][0]*b.x+a[2][1]*b.y+a[2][2]*b.z;
    return  c;
  }

  //Affin変換用
  //VecMtx*VecN　行列の右からベクトルを乗算する
  //行とベクトルの次数があわなければ零で埋まった行列を帰す
  VecN<T>  operator*(VecN<T> &b){
    VecMtx<T> &a = *this;
    int i = a.RowN();  VecN<T> c(i);  if(b.vn != ColumN()) return  c;
    for(int n=0;  n<i;  n++)  {
      c[n]=0; for(int k=0;  k<b.vn; k++)  c[n]+=a[n][k]*b[k];  }
    return  c;
  }

  //VecMtx*VecMtx　行列の乗算　関数版
  //aの列とbの行数があわなければなにもしない
  //c.mul(a,b) は c=a*b と等価
  void  mul(VecMtx<T> &a, VecMtx<T> &b) {
    if(a.Colum != b.Row || a.Row != b.Colum)  return;
    for(int i=0;  i<Row;  i++)  {
      for(int j=0;  j<Colum;  j++)  {
        vp[i][j]=0;
        for(int k=0;  k<Row;  k++)  vp[i][j] += a[i][k] * b[k][j];
      }
    }
  }
  //３行３列の逆行列
  void  Inverse(void) {
    if(Colum==3 && Row==3)  Inverse33();
  }
  void  Inverse33(void)  {
    double  x00 = vp[0][0],   x01 = vp[0][1],   x02 = vp[0][2];
    double  x10 = vp[1][0],   x11 = vp[1][1],   x12 = vp[1][2];
    double  x20 = vp[2][0],   x21 = vp[2][1],   x22 = vp[2][2];
    double  bunbo = x00*(x11*x22-x12*x21)+x01*(x12*x20-x10*x22)+x02*(x10*x21-x11*x20);
    bunbo =1.0/bunbo;
    vp[0][0] = (x11*x22-x12*x21)*bunbo;
    vp[0][1] = (x02*x21-x01*x22)*bunbo;
    vp[0][2] = (x01*x12-x02*x11)*bunbo;
    vp[1][0] = (x12*x20-x10*x22)*bunbo;
    vp[1][1] = (x00*x22-x02*x20)*bunbo;
    vp[1][2] = (x02*x10-x00*x12)*bunbo;
    vp[2][0] = (x10*x21-x11*x20)*bunbo;
    vp[2][1] = (x01*x20-x00*x21)*bunbo;
    vp[2][2] = (x00*x11-x01*x10)*bunbo;
  }

  //転置行列への変換
  void  Transpose(void) {
    if(Colum!=Row)  return; //行と列があわなければなにもしない
    int n,m;
    T c;
    for(n=0;  n<Row;  n++)  {
      for(m=n+1;  m<Colum;  m++)  {
        c = vp[n][m]; vp[n][m] = vp[m][n];  vp[m][n] = c;
      }
    }
  }

  VecMtx<T> & Assign(const  VecMtx<T> &vm)  {
    int r,  c;
    for(r=0;  r<Row && r<vm.Row; r++) {
      for(c=0;  c<Colum && c<vm.Colum; c++) {
        vp[r][c] = vm.vp[r][c];
      }
    }
    return  *this;
  }

  friend  ostream &operator<<(ostream &s, VecMtx<T> &p) {
    for(int n=0;  n<p.Row; )  s << p.vp[n++] << endl;  return  s;
  }

  friend  istream &operator>>(istream &s, VecMtx<T> &p) {
    for(int n=0;  n<p.vn; ) s >> p.vp[n++];
    return  s;
  }

  //以下の二つの関数は基本レベルの変数タイプだけを扱う char int double
  int WriteBin(FILE *fp) {
    int m, k;
    for(int n=0;  n<Row;  n++)  {
      if((k=vp[n].WriteBin(fp))==0)  return  0;  else  m+=k;
    }
    return  m;
  }
  int WriteBin(FILE *fp,  VecMtx<T> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  VecMtx<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {
    int m, k;
    for(int n=0;  n<Row;  n++)  {
      if((k=vp[n].ReadBin(fp))==0)  return  0;  else  m+=k;
    }
    return  m;
  }
  int ReadBin(FILE *fp,  VecMtx<T> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  VecMtx<T> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }
};
//======================================
//--------------------------------------
//row行 col列　直方行列 VecNs を使った行列
template  <class T, int Row, int Colum>
class  VecMs
{
  private:
  public:
  VecNs<T, Colum> v[Row];

  VecMs()  {}

  int RowN(void) { return Row;  }     //行数
  int ColumN(void) { return Colum;  } //列数

  void  Identity(void)  { //単位行列の生成 Row＞Columならエラーを起こす
    for(int n=0;  n<Row;  n++)  { for(int m=0;  m<Colum;  m++)  v[n][m]=(T)0; }
    for(int n=0;  n<Row;  n++)  v[n][n]=1;
  }

  //行列式の値、２か３の正方行列　それ以外の場合は最大値
  double Det2(void) {
    return  v[0][0]*v[1][1]-v[0][1]*v[1][0];
  }

  double Det3(void) {
    return  v[0][0]*v[1][1]*v[2][2]
            +v[0][1]*v[1][2]*v[2][0]
            +v[0][2]*v[1][0]*v[2][1]
            -v[0][2]*v[1][1]*v[2][0]
            -v[0][1]*v[1][0]*v[2][2]
            -v[0][0]*v[1][2]*v[2][1];

  }

  double Det4(void) {
    return
     v[0][0]*(v[1][1]*(v[2][2]*v[3][3]-v[2][3]*v[3][2])
              +v[1][2]*(v[2][3]*v[3][1]-v[2][1]*v[3][3])
              +v[1][3]*(v[2][1]*v[3][2]-v[2][2]*v[3][1]))
    -v[0][1]*(v[1][0]*(v[2][2]*v[3][3]-v[2][3]*v[3][2])
              +v[1][2]*(v[2][3]*v[3][0]-v[2][0]*v[3][3])
              +v[1][3]*(v[2][0]*v[3][2]-v[2][2]*v[3][0]))
    +v[0][2]*(v[1][0]*(v[2][1]*v[3][3]-v[2][3]*v[3][1])
              +v[1][1]*(v[2][3]*v[3][0]-v[2][0]*v[3][3])
              +v[1][3]*(v[2][0]*v[3][1]-v[2][1]*v[3][0]))
    -v[0][3]*(v[1][0]*(v[2][1]*v[3][2]-v[2][2]*v[3][1])
              +v[1][1]*(v[2][2]*v[3][0]-v[2][0]*v[3][2])
              +v[1][2]*(v[2][0]*v[3][1]-v[2][1]*v[3][0]));
  }

  double Det(void) {
    if(Colum==2 && Row==2)  return  Det2();
    else  if(Colum==3 && Row==3)  return  Det3();
    else  if(Colum==4 && Row==4)  return  Det4();
    else  return  0.0;
  }

  //２次元回転行列
  void  Rotate(double p)  {
    v[0][0] = v[1][1] = cos(p);
    v[1][0] = -(v[0][1] = sin(p));
  }

  //任意の回転ベクトルvを軸にする(p)回転行列要素の生成
  //日刊工業新聞　「入門　ＣＡＤ図形計算手法」ｐ３２
  //３行３列以外では最初にIdentityで単位行列にしてから実行する
  void  RotateV(Vec3<T> &v3, double p)  {
    double  c, c1, s;
    Vec3<T> vt = v3;
    c1 = 1.0-(c = cos(p));  s = sin(p);
    if(vt.power()==0.0)  {  Identity(); return; }
    vt.Standard();
    v[0][0]=c1*vt.x*vt.x+c;      v[0][1]=c1*vt.x*vt.y-vt.z*s; v[0][2]=c1*vt.x*vt.z+vt.y*s;
    v[1][0]=c1*vt.x*vt.y+vt.z*s; v[1][1]=c1*vt.y*vt.y+c;      v[1][2]=c1*vt.y*vt.z-vt.x*s;
    v[2][0]=c1*vt.x*vt.z-vt.y*s; v[2][1]=c1*vt.y*vt.z+vt.x*s; v[2][2]=c1*vt.z*vt.z+c;
  }

  void  RotateVA(Vec3<T> &v3)  {  RotateV(v3, abs(v3));  }

  //回転表現 X,Y,Z　要素だけを変化させる　
  //３行３列以外では最初にIdentityで単位行列にしてから実行する
  void  RotateZ(double  x,  double  y)  {
    if(Row != 2 || Colum != 2)  Identity();
    v[0][0] = v[1][1] = x;    v[1][0] = -(v[0][1]= y);
  }

  void  RotateZ(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    v[0][0] = v[1][1] = cos(p);    v[1][0] = -(v[0][1]= sin(p));
  }
  void  RotateX(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    v[1][1] = v[2][2] = cos(p);    v[1][2] = -(v[2][1]= sin(p));
  }
  void  RotateY(double p)  {
    if(Row != 2 || Colum != 2)  Identity();
    v[0][0] = v[2][2] = cos(p);    v[0][2] = -(v[2][0]= sin(p));
  }

  void  SolveC1(void)  // 対角を１にして、係数を第一列に吐き出す
  { int nC1, nC2;  T  tmpCC;
    for(nC1 = 1;  nC1 < Row; nC1++) {
      tmpCC = v[nC1][nC1];
      for(nC2 = 1;  nC2 < Row;  nC2++)  {
        if(nC1 != nC2)  v[nC2] -=  v[nC1]*(v[nC2][nC1]/tmpCC);
      }
    }

    for(nC1=1;  nC1 < Row;  nC1++)  { // 行列の対角要素を１に規格化する
      tmpCC = 1.0/v[nC1][nC1];  v[nC1] *=  tmpCC;
    }
  }

  VecNs<T, Row>& operator[](int n)  { return  v[n];  }

  VecMs<T, Row, Colum> operator*(T d)  {  //定数培
    VecMs<T, Row, Colum> c;
    for(int m=0;  m<Row;  m++)  {    c[m] = v[m]*d;  }  return c;
  }

  Vec2<T> operator*(Vec2<T> &b)
  {
    VecMs<T, Row, Colum> &a = *this;
    Vec2<T> c;
    c.x = a[0][0]*b.x+a[0][1]*b.y;
    c.y = a[1][0]*b.x+a[1][1]*b.y;
    return  c;
  }

  Vec3<T> operator*(Vec3<T> &b)
  {
    VecMs<T, Row, Colum> &a = *this;
    Vec3<T> c;
    c.x = a[0][0]*b.x+a[0][1]*b.y+a[0][2]*b.z;
    c.y = a[1][0]*b.x+a[1][1]*b.y+a[1][2]*b.z;
    c.z = a[2][0]*b.x+a[2][1]*b.y+a[2][2]*b.z;
    return  c;
  }

  //Affin変換用
  //VecMs*VecNs　行列の右からベクトルを乗算する
  //行とベクトルの次数があわなければ零で埋まった行列を帰す
  VecNs<T, Row>  operator*(VecNs<T, Colum> &b){
    VecNs<T, Row> c;
    for(int n=0;  n<RowN();  n++) c[n]=v[n]*b;
    return  c;
  }

  //VecMs*VecMs　行列の乗算　関数版
  //aの列とbの行数があわなければなにもしない
  //c.mul(a,b) は c=a*b と等価
  void  mul(VecMs<T, Row, Colum> &a, VecMs<T, Row, Colum> &b) {
    if(a.ColumN() != b.RowN() || a.RowN() != b.ColumN())  return;
    for(int i=0;  i<RowN();  i++)  {
      for(int j=0;  j<ColumN();  j++)  {
        v[i][j]=0;
        for(int k=0;  k<Row;  k++)  v[i][j] += a[i][k] * b[k][j];
      }
    }
  }
  //VecMs*VecMs　行列の乗算
  //行の次数があわなければ零で埋まった行列を帰す、内部で関数版を使用している
  //代数学と幾何学（学術図書出版)P11
  VecMs<T, Row, Colum> operator*(VecMs<T, Row, Colum> &b)
  { VecMs<T, Row, Colum> &a = *this;
    VecMs<T, Row, Colum> c;    c.mul(a,b);  return c;}

  //３行３列の逆行列
  void  Inverse(void) {
    if(Colum==3 && Row==3)  Inverse33();
  }
  void  Inverse33(void)  {
    double  x00 = v[0][0],   x01 = v[0][1],   x02 = v[0][2];
    double  x10 = v[1][0],   x11 = v[1][1],   x12 = v[1][2];
    double  x20 = v[2][0],   x21 = v[2][1],   x22 = v[2][2];
    double  bunbo = x00*(x11*x22-x12*x21)+x01*(x12*x20-x10*x22)+x02*(x10*x21-x11*x20);
    bunbo =1.0/bunbo;
    v[0][0] = (x11*x22-x12*x21)*bunbo;
    v[0][1] = (x02*x21-x01*x22)*bunbo;
    v[0][2] = (x01*x12-x02*x11)*bunbo;
    v[1][0] = (x12*x20-x10*x22)*bunbo;
    v[1][1] = (x00*x22-x02*x20)*bunbo;
    v[1][2] = (x02*x10-x00*x12)*bunbo;
    v[2][0] = (x10*x21-x11*x20)*bunbo;
    v[2][1] = (x01*x20-x00*x21)*bunbo;
    v[2][2] = (x00*x11-x01*x10)*bunbo;
  }

  //転置行列への変換
  void  Transpose(void) {
    if(ColumN()!=RowN())  return; //行と列があわなければなにもしない
    int n,m;
    T c;
    for(n=0;  n<Row;  n++)  {
      for(m=n+1;  m<Colum;  m++)  {
        c = v[n][m]; v[n][m] = v[m][n];  v[m][n] = c;
      }
    }
  }

  VecMs<T, Row, Colum> &  Assign(const  VecMs<T, Row, Colum> &vm)  {
    for(int r=0;  r<Row; r++) v[r] = vm.v[r];
    return  *this;
  }

  friend  ostream &operator<<(ostream &s, VecMs<T, Row, Colum> &p)  {
    for(int n=0;  n<p.Row; )  s << p.v[n++] << endl;  return  s;
  }

  friend  istream &operator>>(istream &s, VecMs<T, Row, Colum> &p)  {
    for(int n=0;  n<p.Row; ) s >> p.v[n++];
    return  s;
  }


  //以下の二つの関数は基本レベルの変数タイプだけを扱う char int double
  int WriteBin(FILE *fp) {
    int m, k;
    for(int n=0;  n<Row;  n++)  {
      if((k=v[n].WriteBin(fp))==0)  return  0;  else  m+=k;
    }
    return  m;
  }
  int WriteBin(FILE *fp,  VecMs<T, Row, Colum> *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  VecMs<T, Row, Colum> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {
    int m, k;
    for(int n=0;  n<Row;  n++)  {
      if((k=v[n].ReadBin(fp))==0)  return  0;  else  m+=k;
    }
    return  m;
  }
  int ReadBin(FILE *fp,  VecMs<T, Row, Colum> *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  VecMs<T, Row, Colum> *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }
};
//--------------------------------------
template  <class T, int Row, int Colum>
int ReadBin(char *na,  VecMs<T, Row, Colum> *v, int n) {
  return  v->ReadBin(na, v, n);
}
//--------------------------------------
template  <class T, int Row, int Colum>
int WriteBin(char *na, VecMs<T, Row, Colum> *v, int n) {
  return  v->WriteBin(na, v, n);
}
//======================================
//球座標　内部で正規化されるので必ずラジアンで扱うこと
//2000/11/20 修正　x=y=0のときにエラーを起こしてしまう

class   Spherical
{
  public:
  double  r;    //動径
  double  phi;  //XY面内、+Z軸から見てX軸基準に左回りの角度
  double  theta;//XY面から+Z軸への傾き角度

  //code
  Spherical() { r = phi = theta = 0.0;  }
  Spherical(double  rr, double pp, double  tt) {
    r = rr;   phi = pp;  theta = tt;
    adj();
  }

  explicit Spherical(const Vec3<double> &v3) { Assign(v3); }

  //演算  2000/02/03 Hal.T
  //４次であっても３次として計算する
  #ifndef StrictEqual //{
  Spherical& operator=(VecN<double> &src)  {  return  Assign(src);  }
  Spherical& operator=(VecNs<double,4> &src)  {  return  Assign(src);  }
  #endif  //}

  //2000/11/20 修正　x=y=0のときにエラーを起こしてしまう
  Spherical&  Assign(const Vec3<double> &v3) {
    if((r=v3.abs())==0) {
      phi=theta=0;
    } else  {
      phi = (v3.x==0 && v3.y==0)? 0 : atan2(v3.y, v3.x);
      theta = atan2(v3.z, sqrt(v3.x*v3.x+v3.y*v3.y));
    }
    return  *this;
  }

  Spherical&  Assign(VecN<double> &src) {
    Vec3<double> v3;    v3.Assign(src);
    r = v3.abs();       phi = atan2(v3.y, v3.x);
    theta = atan2(v3.z, sqrt(v3.x*v3.x+v3.y*v3.y));
    return  *this;
  }

  Spherical&  Assign(VecNs<double,4> &src)  {
    Vec3<double> v3;    v3.Assign(src);
    r = v3.abs();       phi = atan2(v3.y, v3.x);
    theta = atan2(v3.z, sqrt(v3.x*v3.x+v3.y*v3.y));
    return *this;
  }

  //  r>=0.0, theta =±π/2,phi=０から２πの範囲に収める
  void  adj(void)  {
    int n;
    if(r<0.0) { r=-r; phi+=PI1;  theta=-theta;  } //r>=0.0にする

    if(PID2<theta || theta<-PID2) {   //-π/2<theta<π/2
      //0-2PIの範囲へ
      n = theta/PI2;  if(theta<0.0) n--;
      theta -= n*PI2;
      if(theta>PID2 && theta<=(PI1)) { //90-180
        theta = PI1-theta; phi+=PI1;
      } else  if(theta>PI1 && theta<(PI1+PID2)) { //180-270
        theta = theta-PI1; phi+=PI1;
      } else  if(theta>=(PI1+PID2))  { //270-360
        theta -= PI2;
      }
    }
    if(0<=phi && phi<PI2) return;       //0<phi<2π
    n = phi/PI2;    if(phi<0.0) n--;
    phi -= n*PI2;
  }

  int WriteBin(FILE *fp) {  return  fwrite(&r,  sizeof(double), 3, fp);  }
  int WriteBin(FILE *fp,  Spherical *v, int n) {
    int i=0,  m=0;
    do  { m+=(i=v->WriteBin(fp)); v++;  n--;  } while(n>0 && i>0);
    return  m;
  }
  int WriteBin(char *name,  Spherical *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "wb"))>0) m = WriteBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }

  int ReadBin(FILE *fp) {  return  fread(&r,  sizeof(double), 3, fp);  }
  int ReadBin(FILE *fp,  Spherical *v, int n) {
    int i=0,  m=0;
    do  {
      m+=(i=v->ReadBin(fp)); v++;  n--;
    } while(n>0 && i>0);
    return  m;
  }
  int ReadBin(char *name,  Spherical *v, int n) {
    int m=0;  FILE  *fp;
    if((fp = fopen(name, "rb"))>0) m = ReadBin(fp,  v,  n);
    fclose(fp);
    return  m;
  }
};

//STLのuniqueを使うために関数テンプレートに変えた 2001/2/21 Hal.T 
bool operator==(const Spherical &p, const Spherical &q)  {
  return p.r==q.r && p.phi==q.phi && p.theta==q.theta; }
bool operator!=(const Spherical &p, const Spherical &q)  {
  return p.r!=q.r || p.phi!=q.phi || p.theta!=q.theta; }


//======================================
//数理構造
//球座標＋別のデータ構造
template  <class T>
class   SphericalVector : public  Spherical
{
  public:
  T v;
  //code
  SphericalVector():Spherical() {}

  SphericalVector(T vv, double  rr=0, double pp=0.0, double  tt = 0.0) :
  Spherical(rr, pp, tt), v(vv)  {}
//  Spherical(rr, pp, tt) { v.Assign(vv);  }

  SphericalVector&  Assign(SphericalVector<T> &p) {
    r=p.r;  phi=p.phi;  theta=p.theta;    v=p.v;
    return  *this;
  }

};
//======================================
//１度は実体の寸前まで作成しておけば　配列の時にエラーが出ない
class DumyVector  {
  Vec2<int>   Vec2_int;
  Vec3<int>   Vec3_int;
  VecN<int>   VecN_int;
  VecMtx<int> VecMtx_int;

  Vec2<double>  Vec2_double;
  Vec3<double>  Vec3_double;
  VecN<double>  VecN_double;
  VecNs<double, 4>  VecNs4_double;
  VecMtx<double> VecMtx_double;

};
//======================================
typedef Vec2<int>       V2I;
typedef Vec2<float>     V2F;
typedef Vec2<double>    V2D;
typedef Vec3<int>       V3I;
typedef Vec3<float>     V3F;
typedef Vec3<double>    V3D;
typedef VecNs<double,3> V3Ds;
typedef VecNs<double,4> V4Ds;
typedef VecMtx<int>     VMI;
typedef VecMtx<float>   VMF;
typedef VecMtx<double>  VMD;
/*
VecNs, VecMsの注意書き
一般に静的に確保されるものは宣言時にサイズが確定していなければならない
*/
//======================================

#endif  // __VECTORN_H