#pragma once
#include "math.h"
#include <list>
#include <vector>
#include <algorithm>
extern const double pi, rdpdg;
class plnpos;

//Definitions of projections are in the later half of this file.


//Basic classes
class spcpos {//3-dimensional position
public: double x, y, z;
	//constructors
	spcpos() : x(0), y(0), z(0) {}
	spcpos(double d, double e, double f) : x(d), y(e), z(f) {}
	//operators
	spcpos operator+=(const spcpos &d) { x += d.x; y += d.y; z += d.z; return *this; }
	spcpos operator-=(const spcpos &d) { x -= d.x; y -= d.y; z -= d.z; return *this; }
	spcpos operator*=(const double &d) { x *= d; y *= d; z *= d; return *this; }
	spcpos operator/=(const double &d) { x /= d; y /= d; z /= d; return *this; }
	spcpos operator-() const { return spcpos(-x, -y, -z); }
	spcpos operator+(const spcpos &d) const { return spcpos(x + d.x, y + d.y, z + d.z); }
	spcpos operator-(const spcpos &d) const { return spcpos(x - d.x, y - d.y, z - d.z); }
	spcpos operator*(const double &d) const { return spcpos(x * d, y * d, z * d); }
	spcpos operator/(const double &d) const { return spcpos(x / d, y / d, z / d); }
	//function
	//coordinate for cutting
	plnpos azeqdstX(double f) const;
	//XYZ <-> spherical coordinate, (0 degree N,0 degree E)=X-axis, N-pole = Z-axis
	spcpos XYsphpos() const { return spcpos(cos(y) * cos(x), cos(y) * sin(x), sin(y)) * z; }
	spcpos XYsphcrd() const {
		double r = sqrt(x * x + y * y + z * z), s = sqrt(x * x + y * y);
		if (s == 0) { return spcpos(0, atan2(z, s), r); } else { return spcpos(atan2(y, x), atan2(z, s), r); }
	}
	//XYZ <-> cylindrical coordinate, (0 degree N,0 degree E)=X-axis, N-pole = Z-axis
	spcpos XYcylpos() const { return spcpos(y * cos(x), y * sin(x), z); }
	spcpos XYcylcrd() const {//xyz to cylindrical coordinate
		double s = sqrt(x * x + y * y);
		if (s == 0) { return spcpos(0, s, z); } else { return spcpos(atan2(y, x), s, z); }
	}
};

class plnpos {//s-dimensional position
public: double x, y;
	//constructors
	plnpos() : x(0), y(0) {}
	plnpos(double d, double e) : x(d), y(e) {}
	//operators
	plnpos operator+=(const plnpos &d) { x += d.x; y += d.y; return *this; }
	plnpos operator-=(const plnpos &d) { x -= d.x; y -= d.y; return *this; }
	plnpos operator*=(const double &d) { x *= d; y *= d; return *this; }
	plnpos operator/=(const double &d) { x /= d; y /= d; return *this; }
	plnpos operator-() const { return plnpos(-x, -y); }
	plnpos operator+(const plnpos &d) const { return plnpos(x + d.x, y + d.y); }
	plnpos operator-(const plnpos &d) const { return plnpos(x - d.x, y - d.y); }
	plnpos operator*(const double &d) const { return plnpos(x * d, y * d); }
	plnpos operator/(const double &d) const { return plnpos(x / d, y / d); }
	//functions
	double sqlength() const { return  x * x + y * y; }
	double inproduct(plnpos &d) const { return x * d.x + y * d.y; }
	double vcproduct(plnpos &d) const { return x * d.y - y * d.x; }
	plnpos dsply(double f, plnpos s) { return plnpos(x, -y) * f / 2.0 + s; }//to convert to display coordinate
	//coordinate for cutting
	spcpos RazeqdstX(double f) const;
	plnpos polar2XY() { return plnpos(cos(y), sin(y)) * x; }
	plnpos XY2polar() {
		double r = sqrt(this->sqlength());
		if (r == 0) { return plnpos(0.0, 0.0); }
		else { return plnpos(r, atan2(y, x)); }
	}
	//operators as complex numbers
	plnpos operator*(plnpos &d) { return plnpos(x * d.x - y * d.y, y * d.x + x * d.y); }
	plnpos operator/(plnpos &d) { return plnpos(x * d.x + y * d.y, y * d.x - x * d.y) / d.sqlength(); }
	plnpos complexsqrt() {
		double r = sqrt(sqlength());
		plnpos z = plnpos(sqrt((r + x)/ 2.0), sqrt((r - x) / 2.0));
		if (y < 0) { z.y *= -1; }
		return z;
	}
	//the angle between 2 vectors
	double angle(plnpos v1) {
		double cs = this->inproduct(v1), sn = this->vcproduct(v1);
		if (cs == 0 && sn == 0) { return 0; } else { return atan2(sn, cs); }
	}
	//double dsttoline(plnpos v1, plnpos v2) {
	//	plnpos d0 = v2 - v1, d1 = *this - v1, d2 = *this - v2;
	//	double r0 = d0.sqlength();
	//	double t = d1.inproduct(d0) / r0;
	//	if (t > 1) { return sqrt(d2.sqlength()); }
	//	else if (t < 0) { return sqrt(d1.sqlength()); }
	//	else { return abs(d2.vcproduct(d1)) / sqrt(r0); }
	//}
};

class mtrx3 {//matrix of order 3 for 3-dimensional rotations
public: spcpos v0, v1, v2;
	//constructor
	mtrx3() {}
	mtrx3(spcpos d0, spcpos d1, spcpos d2) : v0(d0), v1(d1), v2(d2) {}
	mtrx3(spcpos d) {
		//(0 degree N, 0 degree E)=X-axis, (0 degree N, 90 degree E)=Y-axis, N-pole=Z-axis
		//then this matrix rotates a position(lng, lat, rot) to X-axis.
		//This means lng-clockwise rotation around Z-axis, lat-clockwise rotation around Y-axis, and rot-clockwise rotation around X-axis
		double csx = cos(d.x), snx = sin(d.x), csy = cos(d.y), sny = sin(d.y), csz = cos(d.z), snz = sin(d.z);
		v0 = spcpos(csx * csy, -snx * csz - csx * sny * snz,  snx * snz - csx * sny * csz);
		v1 = spcpos(snx * csy,  csx * csz - snx * sny * snz, -csx * snz - snx * sny * csz);
		v2 = spcpos(sny, csy * snz, csy * csz);
	}
	//operator
	spcpos operator*(const spcpos &d) const { return v0 * d.x + v1 * d.y + v2 * d.z; }
	mtrx3 operator*(const mtrx3 &d) const { return mtrx3(*this * d.v0, *this * d.v1, *this * d.v2); }
	mtrx3 operator*(const double &d) const { return mtrx3(v0 * d, v1 * d, v2 * d); }
	//function
	mtrx3 transpose() {
		spcpos u0 = spcpos(v0.x, v1.x, v2.x), u1 = spcpos(v0.y, v1.y, v2.y), u2 = spcpos(v0.z, v1.z, v2.z);
		return mtrx3(u0, u1, u2);
	}
	//const
	static const mtrx3 X2X, N2X, S2X, W2X, E2X, A2X, ATX, R045, R090, R270, R315;
};

//2nd level
class crosspnt {// for ordering cut points around origin
public: int name, next, pntnum, pgnum, edge, drct, ch;
	double pnt;//parameter on the edge
	plnpos sstn, plcrd; spcpos ssph, *sph;
	crosspnt(int n, int x, int pg, int a, double b, int d, plnpos p, spcpos s) : name(n), next(x), pntnum(0), pgnum(pg), edge(a), pnt(b), drct(d), ch(0), ssph(s) {
		sstn = p; plcrd = p.XY2polar(); sph = nullptr;
	}
};

class plygn {
public: int type, pntnum, pntname;//type: 1:polygon, 2:lake, 3:line, 4:indicatrix, 5: hole of indicatrix, 10:sea, 11:hole of sea
	int incl, inclname; spcpos *sph;
	//constructor
	plygn(int t, int n, spcpos *p1, int d, int e, int f) : type(t), pntnum(n), sph(p1)  {
		pntname = d; incl = e; inclname = f;
	}
	//function
	//counting the number of rotation around the origin.
	double rotation(mtrx3 mt1, double f, double xcut) {
		double aa = 0; spcpos *tmpsph = sph, sv0, sv1 = mt1 * *tmpsph; plnpos v0, v1 = sv1.azeqdstX(f);
		for (int i = 0; i < pntnum; ++i) {
			sv0 = sv1; v0 = v1;
			if (i == pntnum - 1) { tmpsph = sph; } else { ++tmpsph; }
			sv1 = mt1 * *tmpsph; v1 = sv1.azeqdstX(f);
			if (sv0.x > xcut && sv1.x > xcut) { aa += v0.angle(v1); }
		}
		//if (incl % 2 == 1) { aa *= -1.0; }
		return aa;
	}
	//size of area
	double sizeofarea(mtrx3 mt1, double f, double xcut) {
		double ss = 0; spcpos *tmpsph = sph, sv0, sv1 = mt1 * *tmpsph; plnpos v0, v1 = sv1.azeqdstX(f);
		for (int i = 0; i < pntnum; ++i) {
			sv0 = sv1; v0 = v1;
			if (i == pntnum - 1) { tmpsph = sph;  } else { ++tmpsph; }
			sv1 = mt1 * *tmpsph; v1 = sv1.azeqdstX(f);
			if (sv0.x > xcut && sv1.x > xcut) { ss += v0.vcproduct(v1); }
		}
		//if (incl % 2 == 1) { ss *= -1.0; }
		return ss;
	}
	void linecutter(mtrx3 mt1, mtrx3 mt2, spcpos cut, double xcut, std::list<plygn> &line, int &linemax) {
		int cpre = -1, cpre2 = -1, cpst = -1, cmin = -1, cmin2 = -1, cmax = -1;
		double r0, r1, r2, t0, t1;
		double sqcutting = cut.x * cut.x;
		spcpos *newsph, *tmpsph = sph, vpre, vpre2, vpst, sv0, sv1 = mt1 * *tmpsph;
		plnpos v0, v1 = sv1.azeqdstX(cut.y), vd;
		for (int i = 0; i < pntnum; ++i) {
			//check conditions of cutting
			if (i < pntnum - 1) {
				sv0 = sv1; v0 = v1; ++tmpsph;
				sv1 = mt1 * *tmpsph; v1 = sv1.azeqdstX(cut.y);
				if (sv0.x > xcut || sv1.x > xcut) {
					vd = v1 - v0; r2 = vd.sqlength();
					r0 = -vd.inproduct(v0) / r2;
					r1 = v1.vcproduct(v0) / r2; r1 = sqcutting / r2 - r1 * r1;
					if (r1 > 0) {//the end point
						t0 = r0 - sqrt(r1); t1 = r0 + r0 - t0;
						if (t0 > 0.0 && t0 < 1.0) {
							vpst = mt2 * (vd * t0 + v0).RazeqdstX(cut.y); cmax = i; cpst = 1;
							if (cmin < 0) { cmin = i; cpre = 0; }
						} else if (t1 <= 0.0 || t0 >= 1.0) {
							if (cmin < 0) { cmin = i; cpre = 0; }
						} else if (cmin >= 0) { cmax = i; cpst = 0; }
						//preparation for the next begining point
						if (t1 > 0.0 && t1 < 1.0) {
							vpre2 = mt2 * (vd * t1 + v0).RazeqdstX(cut.y); cmin2 = i + 1; cpre2 = 1;
						}
					} else if (cmin < 0) { cmin = i; cpre = 0; }
				} else if (cmin < 0) { cmin = i; cpre = 0; }
			} else { cmax = i; cpst = 0; }//the end of the loop
		 //output new line
			if (cmin >= 0 && cmax >= 0) {
				newsph = new spcpos[cmax - cmin + cpre + cpst + 1];
				if (cpre > 0) {
					newsph[0] = vpre;
				}
				if (cpst > 0) {
					newsph[cmax - cmin + cpre + cpst] = vpst;
				}
				for (int j = 0; j <= cmax - cmin; ++j) {
					newsph[j + cpre] = sph[j + cmin];
				}
				++linemax;
				line.push_back(plygn(3, cmax - cmin + cpre + cpst + 1, newsph, linemax, 0, 0));
				cmin = -1; cmax = -1; cpre = -1; cpst = -1;
			}
			if (cmin2 >= 0) {//the begining point
				vpre = vpre2; cmin = cmin2; cpre = cpre2;
				cmin2 = -1; cpre2 = -1;
			}
		}
	}
};

class plygnset {
public: std::list<plygn> plygon; int pntname, type, plych;//plych:2:no polygon but all of sphere, 1:polygons exist 0:no polygon & no part
	//constructor
	plygnset(plygn p, int n, int t) { plygon.push_back(p); pntname = n; type = t; plych = 1; }
	plygnset(int n, int t, int c) { pntname = n; type = t; plych = c; }
	//functions
	int rotation(mtrx3 mt1, double f, double xcut) {
		int rot = 0; double r;
		for (auto itr = plygon.begin(); itr != plygon.end(); itr++) {
			r = itr->rotation(mt1, f, xcut);
			if (r > 2.0) { ++rot; } else if (r < -2.0) { --rot; }
		}
		return rot;
	}
	double sizeofarea(mtrx3 mt1, double f, double xcut) {
		double ss = 0;
		for (auto itr = plygon.begin(); itr != plygon.end(); itr++) { ss += itr->sizeofarea(mt1, f, xcut); }
		return ss;
	}
	void plygncutter(mtrx3 mt1, mtrx3 mt2, double sh1, spcpos cut, double xcut) {
		double sqcutting = cut.x * cut.x, xcut2 = -cos(cut.z * 0.5);
		std::vector<crosspnt> crs1;
		plnpos v0, v1, vd; spcpos sv0, sv1, svp, *tmpsph;
		double r0, r1, r2, t0, t1;

		//cnt:number of cutting point
		int cnt = 0, pcnt, pnum = plygon.size();
		for (int j = 0; j < pnum; ++j) {
			pcnt = cnt;
			//ch:check for cutting exactly begining or ending point, eh:check for completely containing(outsie:1,inside:-1)
			int ch = 0, eh = -1;
			tmpsph = plygon.begin()->sph;
			sv1 = mt1 * *(tmpsph + plygon.begin()->pntnum - 1); v1 = sv1.azeqdstX(cut.y);
			for (int i = -1; i < plygon.begin()->pntnum; ++i) {//calculating cross point, in:+ , out:-
				sv0 = sv1; v0 = v1;
				if (i == plygon.begin()->pntnum - 1) { sv1 = mt1 * *(plygon.begin()->sph);  v1 = sv1.azeqdstX(cut.y); }
				else { sv1 = mt1 * *tmpsph; v1 = sv1.azeqdstX(cut.y); ++tmpsph; }
				if (sv0.x > xcut || sv1.x > xcut) {
					vd = v1 - v0; r2 = vd.sqlength();
					r0 = -vd.inproduct(v0) / r2;
					r1 = v1.vcproduct(v0) / r2; r1 = sqcutting / r2 - r1 * r1;
					if (r1 > 0) {
						t0 = r0 - sqrt(r1); t1 = r0 + r0 - t0;
						if (t0 <= 0.0 && t1 > 0.0) {
							if (ch > 0) {
								svp = mt2 * v0.RazeqdstX(cut.y);
								crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, 5, v0, svp));
								//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, 5, v0, svp));
								++cnt; eh = 0;
							}
						} else {
							if (ch < 0) {
								svp = mt2 * v0.RazeqdstX(cut.y);
								crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
								//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
								++cnt; eh = 0;
							}
							if (t0 > 0.0 && t0 < 1.0 && ch != 0) {
								svp = mt2 * (vd * t0 + v0).RazeqdstX(cut.y);
								crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, t0, 5, vd * t0 + v0, svp));
								//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, t0, 5, vd * t0 + v0, svp));
								++cnt; eh = 0;
							}
							eh = abs(eh);
						}
						if (t0 < 1.0 && t1 >= 1.0) {
							ch = -1;
						} else {
							if (t1 > 0.0 && t1 < 1.0 && ch != 0) {
								svp = mt2 * (vd * t1 + v0).RazeqdstX(cut.y);
								crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, t1, -5, vd * t1 + v0, svp));
								//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, t1, -5, vd * t1 + v0, svp));
								++cnt; eh = 0;
							}
							ch = 1; eh = abs(eh);
						}
					} else {
						if (ch < 0) {
							svp = mt2 * v0.RazeqdstX(cut.y);
							crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
							//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
							++cnt; eh = 0;
						}
						ch = 1; eh = abs(eh);
					}
				} else {
					if (ch < 0) {
						svp = mt2 * v0.RazeqdstX(cut.y);
						crs1.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
						//crs2.push_back(crosspnt(cnt, cnt + 1, plygon.begin()->pntname, i, 0, -5, v0, svp));
						++cnt; eh = 0;
					}
					ch = 1; eh = abs(eh);
				}
			}
			if (eh == -1) { delete plygon.begin()->sph; } 
			else if (eh == 1) { plygon.push_back(*(plygon.begin())); }
			else {
				crs1[cnt - 1].next = pcnt;//closing a loop
				int i1, num, cmin, cmax;
				for (int i0 = pcnt; i0 < cnt; ++i0) {//remaining if out
					i1 = (i0 == cnt - 1 ? pcnt : i0 + 1);
					if (crs1[i0].drct < 0) {
						cmin = crs1[i0].edge; cmax = crs1[i1].edge + (crs1[i1].pnt == 0.0 ? 0 : 1);
						if (cmin >= cmax) { cmax += plygon.begin()->pntnum; }
						num = cmax - cmin;
						crs1[i0].pntnum = num;
						crs1[i0].sph = new spcpos[num]; tmpsph = crs1[i0].sph;
						*tmpsph = crs1[i0].ssph; ++tmpsph;
						for (int j = cmin + 1; j < cmax; ++j) {
							*tmpsph = plygon.begin()->sph[j % plygon.begin()->pntnum]; ++tmpsph;
						}
					}
				}
				delete plygon.begin()->sph; 
			}
			plygon.pop_front();
		}
		
		if (cnt > 0) {//if cross points exist
			if (sh1 > 0) {
				std::sort(crs1.begin(), crs1.end(), [](crosspnt c, crosspnt d) {return c.plcrd.y < d.plcrd.y; });//sorting by azimuth
				if (crs1[0].drct < 0) { crs1[0].plcrd.y += pi * 2.0; }//adjust if the first point is "out"
			} else {
				std::sort(crs1.begin(), crs1.end(), [](crosspnt c, crosspnt d) {return c.plcrd.y > d.plcrd.y; });
				if (crs1[0].drct < 0) { crs1[0].plcrd.y -= pi * 2.0; }
			}
			int i1, num; double t;
			for (int i0 = 0; i0 < cnt; ++i0) {//truncating if in
				i1 = (i0 == cnt - 1 ? 0 : i0 + 1);
				if (crs1[i0].drct > 0) {
					crs1[i0].next = crs1[i1].name;
					num = (int)floor((crs1[i1].plcrd.y - crs1[i0].plcrd.y) * sh1 / cut.z) + 1;
					crs1[i0].pntnum = num;
					crs1[i0].sph = new spcpos[num]; tmpsph = crs1[i0].sph;
					t = crs1[i0].plcrd.y;
					for (int j = 0; j < num; ++j) {
						*tmpsph = mt2 * plnpos(cut.x, t).polar2XY().RazeqdstX(cut.y);
						++tmpsph; t += cut.z * sh1;
					}
				}
			}
			std::sort(crs1.begin(), crs1.end(), [](crosspnt c, crosspnt d) {return c.name < d.name; });//resort
			int k;
			spcpos *tmpsph2, *newsph;
			for (int i = 0; i < cnt; ++i) {//connect and output new polygons
				if (crs1[i].ch == 0) {
					num = 0; k = i;
					do { num += crs1[k].pntnum; k = crs1[k].next; } while (k != i);
					newsph = new spcpos[num]; tmpsph = newsph;
					do {
						tmpsph2 = crs1[k].sph;
						for (int j = 0; j < crs1[k].pntnum; ++j) { *tmpsph = *tmpsph2; ++tmpsph; ++tmpsph2; }
						delete[] crs1[k].sph;
						k = crs1[k].next; crs1[k].ch = 1;
					} while (k != i);
					plygon.push_back(plygn(type, num, newsph, pntname, 0, pntname));
				}
			}
		} else {
			pnum = plygon.size();
			if (plych == 1 && pnum == 0) { plych = 0; }
			if (plych > 0) {
				int rt = rotation(mt1, cut.y, xcut2); double ss = sizeofarea(mt1, cut.y, xcut2); 
				int sch = (rt + (ss > 0.0 ? 0 : 1) + (sh1 > 0 ? 0 : 1)) % 2;
				if (plych == 2 || sch == 0) {
					//make a hole
					int num = (int)floor(pi * 2.0 / cut.z) + 1; double t = 0;
					spcpos *newsph = new spcpos[num]; tmpsph = newsph;
					for (int j = 0; j < num; ++j) {
						*tmpsph = mt2 * plnpos(cut.x, t).polar2XY().RazeqdstX(cut.y);
						++tmpsph; t += cut.z * sh1;
					}
					plygon.push_back(plygn(type, num, newsph, pntname, 0, pntname));
					plych = 1;
				}
			}
		}
	}
};


//===============
//projection data
// A definition of a new projection needs two parts.
//  New function.
//  New "case" of the switch statement in the constructor.
//===============
class projection {
public: plnpos(projection::*prjfnc)(spcpos);
		mtrx3 mt1, mt2, mcut1[16], mcut2[16];//1:the position to X-axis, 2:X-axis to the position
		spcpos dflat, dcut[16];//x:raidus, y:factor, z:interval
		int numcut;
		double fctr = 1.0, n = 1.0, c = 1.0;
		static constexpr double exceptionalmargin = 0.000001;
		//constructor
		projection(int prj, spcpos clnglat, spcpos dlat, double lc, double pc, double inv) {
			dflat = dlat; double lctan2 = lc * 2.0 / pi; mt1 = mtrx3(clnglat); mt2 = mt1.transpose();

			//setting a function of projection and matrice to cut the sphere
			//Be careful around points of infinite scale (such as antipodes in finite azimuthal projections, poles in finite cylindrical projections, etc. )
			//If the cutting collapses, try to change the order of matrice.
			switch (prj) {
			case 1:
				prjfnc = &projection::cylindricalequidistant;//Cylindrical projection
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line from N-pole to S-pole
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);//S-pole
				break;
			case 2:
				prjfnc = &projection::cylindricalequalarea;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 3:
				prjfnc = &projection::mercator;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 11:
				prjfnc = &projection::equidistantconic;//Conic projection
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line from N-pole to S-pole
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);//S-pole
				n = (cos(dflat.x) - cos(dflat.y)) / (dflat.y - dflat.x);//factor for azimuth
				c = (dflat.y * cos(dflat.x) - dflat.x * cos(dflat.y)) / (cos(dflat.x) - cos(dflat.y));//factor defined only by dflnglat
				fctr = abs(c) + pi / 2.0;//(factor for size) = (radial coordinate for N-pole or east edge)
				break;
			case 12:
				prjfnc = &projection::albersequalareaconic;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				n = (sin(dflat.x) + sin(dflat.y)) / 2.0;
				c = cos(dflat.x); c = c * c + n * sin(dflat.x) * 2.0;
				fctr = sqrt(c + abs(n) * 2.0) / abs(n);
				break;
			case 13:
				prjfnc = &projection::lambertequalareaconic;
				numcut = 2;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::S2X; dcut[1] = spcpos(pc, 1.0, inv);
				n = (1.0 + sin(dflat.x)) / 2.0;
				fctr = sqrt(2.0 / n);
				break;
			case 14:
				prjfnc = &projection::lambertconformalconic;
				numcut = 2;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				if (dflat.x >= 0) {
					mcut1[1] = mtrx3::S2X; dcut[1] = spcpos(pc, 1.0, inv);
				}
				else {
					mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				}
				n = sin(dflat.x);
				c = cos(dflat.x) * exp(log(tan(pi / 4.0 + dflat.x / 2.0)) * n) / n;
				break;
			case 21:
				prjfnc = &projection::azimuthalequidistant;//Azimuthal projection
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pc, 1.0, inv);//the antipodes
				break;
			case 22:
				prjfnc = &projection::azimuthalequalarea;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pc, 1.0, inv);
				break;
			case 23:
				prjfnc = &projection::orthographic;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, 1.0, inv);//eliminating the antipode hemisphere
				break;
			case 24:
				prjfnc = &projection::gnomonic;
				numcut = 2;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0 + exceptionalmargin, 1.0, inv);//eliminating the antipode hemisphere
				mcut1[1] = mtrx3::A2X; dcut[1] = spcpos(pc, 1.0, inv);//intentional restriction for an infinite projection
				break;
			case 25:
				prjfnc = &projection::stereographic;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pc, 1.0, inv);
				break;
			case 41:
				prjfnc = &projection::eckert1;//Pseudo cylindrical projection with lengthy poles
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);//S-pole
				break;
			case 42:
				prjfnc = &projection::eckert2;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 43:
				prjfnc = &projection::eckert3;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 44:
				prjfnc = &projection::eckert4;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 45:
				prjfnc = &projection::eckert5;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 46:
				prjfnc = &projection::eckert6;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 47:
				prjfnc = &projection::kavrayskiy7;
				numcut = 3;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				break;
			case 61:
				prjfnc = &projection::sinusoidal;//Pseudo cylindrical, conic or azimuthal projection with pointed poles
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line
				break;
			case 62:
				prjfnc = &projection::mollweide;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 63:
				prjfnc = &projection::werner;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 64:
				prjfnc = &projection::bonne;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 65:
				prjfnc = &projection::loximuthal;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 66:
				prjfnc = &projection::americanpolyconic;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 67:
				prjfnc = &projection::rectangularpolyconic;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				if (dflat.x == 0) { fctr = pi; }
				else { fctr = tan(pi * sin(dflat.x) / 2.0) / sin(dflat.x) * 2.0; };
				break;
			case 68:
				prjfnc = &projection::aitoff;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 69:
				prjfnc = &projection::hammer;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 70:
				prjfnc = &projection::lagrange;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 71:
				prjfnc = &projection::eisenlohr;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 72:
				prjfnc = &projection::augustepicycloidal;
				numcut = 1;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				break;
			case 81:
				prjfnc = &projection::interruptedhomolosine;//interrupted projection
				numcut = 7;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);
				mcut1[3] = mtrx3(spcpos(80.0, -45.0, 0) * rdpdg); dcut[3] = spcpos(pi / 4.0, lctan2, inv);
				mcut1[4] = mtrx3(spcpos(-20.0, -45.0, 0) * rdpdg); dcut[4] = spcpos(pi / 4.0, lctan2, inv);
				mcut1[5] = mtrx3(spcpos(-100.0, -45.0, 0) * rdpdg); dcut[5] = spcpos(pi / 4.0, lctan2, inv);
				mcut1[6] = mtrx3(spcpos(-40.0, 45.0, 0) * rdpdg); dcut[6] = spcpos(pi / 4.0, lctan2, inv);
				break;
			case 101:
				prjfnc = &projection::craigretroazimuthal;//front side
				numcut = 4;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[3] = mtrx3(spcpos(0.0, dflat.x, 0.0)) * mtrx3::A2X; dcut[3] = spcpos(pi / 2.0, 1.0, inv);//eliminating the antipode hemisphere of the definiton point
				break;
			case 102:
				prjfnc = &projection::craigretroazimuthal;//back side
				numcut = 4;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//the antipode line
				mcut1[1] = mtrx3::N2X; dcut[1] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[2] = mtrx3::S2X; dcut[2] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[3] = mtrx3(spcpos(0.0, dflat.x, 0.0)); dcut[3] = spcpos(pi / 2.0, 1.0, inv);//eliminating the hemisphere of the definition point

				break;
			case 103:
				prjfnc = &projection::hammerretroazimuthal;//front side
				numcut = 3;
				mcut1[0] = mtrx3::N2X; dcut[0] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[1] = mtrx3::S2X; dcut[1] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[2] = mtrx3::A2X; dcut[2] = spcpos(pi / 2.0, 1.0, inv);//eliminating the antipode hemisphere
				break;
			case 104:
				prjfnc = &projection::hammerretroazimuthal;//back side
				numcut = 4;
				mcut1[0] = mtrx3::N2X; dcut[0] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[1] = mtrx3::S2X; dcut[1] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[2] = mtrx3::X2X; dcut[2] = spcpos(pi / 2.0, 1.0, inv);//eliminating the frontside hemisphere
				mcut1[3] = mtrx3(spcpos(0.0, dflat.x, 0.0)) * mtrx3::A2X; dcut[3] = spcpos(pc, 1.0, inv);//the antipodes of the definiton point
				break;
			case 105:
				prjfnc = &projection::littrow;
				numcut = 3;
				mcut1[0] = mtrx3::N2X; dcut[0] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[1] = mtrx3::S2X; dcut[1] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[2] = mtrx3::A2X; dcut[2] = spcpos(pi / 2.0, 1.0, inv);//eliminating the antipode hemisphere			
				break;
			case 111:
				prjfnc = &projection::peircequincuncial;
				numcut = 6;
				mcut1[0] = mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2, inv);//vertical line of the antipode
				mcut1[1] = mtrx3::ATX; dcut[1] = spcpos(pi / 2.0, lctan2, inv);//horizontal line of the antipode
				mcut1[2] = mtrx3::N2X; dcut[2] = spcpos(pc, 1.0, inv);//N-pole
				mcut1[3] = mtrx3::S2X; dcut[3] = spcpos(pc, 1.0, inv);//S-pole
				mcut1[4] = mtrx3::W2X; dcut[4] = spcpos(pc, 1.0, inv);//W-pole
				mcut1[5] = mtrx3::E2X; dcut[5] = spcpos(pc, 1.0, inv);//E-pole
				break;
			case 112:
				prjfnc = &projection::guyou;
				numcut = 2;
				mcut1[0] = mtrx3::R045 * mtrx3::A2X; dcut[0] = spcpos(pi / 2.0, lctan2 * 20.0, inv);
				mcut1[1] = mtrx3::R045 * mtrx3::ATX; dcut[1] = spcpos(pi / 2.0, lctan2 * 20.0, inv);
				break;
			case 200:
				prjfnc = &projection::test;//test1
				numcut = 2;
				mcut1[0] = mtrx3::R090 * mtrx3::W2X; dcut[0] = spcpos(pi / (2.0 - exceptionalmargin), lctan2, inv);
				mcut1[1] = mtrx3::A2X; dcut[1] = spcpos(pc, 1.0, inv);
				break;
			case 201:
				prjfnc = &projection::azimuthalequidistant;//test2
				numcut = 4;
				mcut1[0] = mtrx3(spcpos(-90.0, 0.0, 2.0) * rdpdg); dcut[0] = spcpos(pi / 2.0, lc * 250.0, inv);
				mcut1[1] = mtrx3(spcpos( 0.0,-87.0, 0.0) * rdpdg); dcut[1] = spcpos(pc * 150.0, 1.0, inv);
				mcut1[2] = mtrx3(spcpos( 0.0, 93.0, 0.0) * rdpdg); dcut[2] = spcpos(pc * 150.0, 1.0, inv);
				mcut1[3] = mtrx3::A2X; dcut[3] = spcpos(pc, 1.0, inv);//the antipodes
				break;
			case 202:
				prjfnc = &projection::azimuthalequalarea;
				numcut = 3;
				mcut1[0] = mtrx3::N2X;  dcut[0] = spcpos(pi / 2.0, lctan2 * 100.0, inv);
				mcut1[1] = mtrx3::X2X;  dcut[1] = spcpos(pc * 100.0, 1.0, inv);
				mcut1[2] = mtrx3::A2X;  dcut[2] = spcpos(pc * 100.0, 1.0, inv);
				break;
			default:
				prjfnc = &projection::cylindricalequidistant;
				numcut = 3;
				mcut1[0] = mtrx3::A2X;  dcut[0] = spcpos(pi / 2.0, lctan2, inv);
				mcut1[1] = mtrx3::N2X;  dcut[1] = spcpos(pc, 1.0, inv);
				mcut1[2] = mtrx3::S2X;  dcut[2] = spcpos(pc, 1.0, inv);
			}
			for (int i = 0; i < numcut; ++i) { mcut2[i] = mcut1[i].transpose(); }
		}
		//inverse function
		double inversefnc(double y, double sx, double(*fnc)(double), double(*dif)(double)) {
			double tx = sx, dy = fnc(tx) - y, dd = dif(tx);
			while (abs(dy) > exceptionalmargin && abs(dd) * 32.0 > exceptionalmargin) { tx -= dy / dd; dy = fnc(tx) - y; dd = dif(tx); }
			return tx;
		}

		//formulae of projections//
		//Input
		//spcpos ipos [x,y,z]:XYZ-coordinate of the sphere.
		//                     X-axis is your side, Y-axis is right-ward, Z-axis is upward.
		//                     The centering point come at X-axis.
		//ipos.XYsphcrd() = (longitude, latitude, radial): Sphere coordinate.
		//                     X-axis is 0 degree N, 0 degree E. Y-axis is 0 degree, 90 degree E, Z-axis is N-pole.
		//ipos.XYcylcrd() = (longitude, cos(latitude), sin(latitude)): Cylindrical coordinate.
		//ipos.z = sin(latitude), so you can write formulae easily on XYZ-coordinate in some projections, especially in azimuthal projections.
		//
		//Output
		//X-axis is right-ward, Y-axis is upward.
		//output coordinates are adjusted to -1.0 <= X <= 1.0 (except infinte projections and some complex projections).
		//The range of Y is varied by projection.

		plnpos cylindricalequidistant(spcpos ipos) { spcpos sph = ipos.XYsphcrd(); return plnpos(sph.x, sph.y) / pi; }
		plnpos cylindricalequalarea(spcpos ipos) { spcpos cyl = ipos.XYcylcrd(); return plnpos(cyl.x, cyl.z) / pi; }
		plnpos mercator(spcpos ipos) { spcpos cyl = ipos.XYcylcrd(); return plnpos(cyl.x, log(sqrt((1.0 + cyl.z) / (1.0 - cyl.z)))) / pi; }
		plnpos equidistantconic(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			double r = c - sph.y, s = n * sph.x;
			return plnpos(r * sin(s), -r * cos(s)) / fctr;
		}
		plnpos albersequalareaconic(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd();
			double r = sqrt(c - n * cyl.z * 2.0) / n, s = n * cyl.x;
			return plnpos(r * sin(s), -r * cos(s)) / fctr;
		}
		plnpos lambertequalareaconic(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd();
			double r = sqrt((1.0 - cyl.z) / n), s = n * cyl.x;
			return plnpos(r * sin(s), -r * cos(s)) / fctr;
		}
		plnpos lambertconformalconic(spcpos ipos) {//with one parameter
			spcpos cyl = ipos.XYcylcrd();
			double r = c / exp(log(sqrt((1.0 + cyl.z) / (1.0 - cyl.z))) * n), s = n * cyl.x;
			return plnpos(r * sin(s), -r * cos(s)) / pi;
		}
		plnpos azimuthalequidistant(spcpos ipos) {
			plnpos opos(ipos.y, ipos.z); double r = sqrt(opos.sqlength());
			if (r == 0) { return plnpos(0, 0); }
			else { return opos * atan2(r, ipos.x) / r / pi; }
		}
		plnpos azimuthalequalarea(spcpos ipos) {
			plnpos opos(ipos.y, ipos.z); double r2 = opos.sqlength();
			if (r2 == 0) { return plnpos(0, 0); }
			else { return opos * sqrt((1.0 - ipos.x) / 2.0 / r2); }
		}
		plnpos orthographic(spcpos ipos) { return plnpos(ipos.y, ipos.z); }
		plnpos gnomonic(spcpos ipos) { return plnpos(ipos.y, ipos.z) / ipos.x; }
		plnpos stereographic(spcpos ipos) { return plnpos(ipos.y, ipos.z) / (ipos.x + 1.0); }
		plnpos sinusoidal(spcpos ipos) { spcpos cyl = ipos.XYcylcrd(); return plnpos(cyl.x * cyl.y, atan2(cyl.z, cyl.y)) / pi; }
		plnpos werner(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd(); double r = pi / 2.0 - sph.y, s = cos(sph.y) / r * sph.x;
			return plnpos(r * sin(s), pi * 15.94 / 45.0 - r * cos(s)) * 1.545 / pi;
		}
		plnpos bonne(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			double r = dflat.x - sph.y + 1.0 / tan(dflat.x), s = cos(sph.y) / r * sph.x;
			return plnpos(r * sin(s), 1.0 / tan(dflat.x) - r * cos(s) + dflat.x) / pi;
		}
		plnpos mollweide(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd(); double y = pi * cyl.z;
			double t = inversefnc(y, cyl.z / 2.5, [](double s) { return sin(s * 2.0) + s * 2.0; }, [](double s) { return cos(s * 2.0) * 2.0 + 2.0; });
			return plnpos(cyl.x * cos(t) / pi, sin(t) / 2.0);
		}
		plnpos homolosine(spcpos ipos) {//not interrupted
			if (abs(ipos.z) > 0.65258) {
				return mollweide(ipos) / 1.1107207 + plnpos(0, (ipos.z > 0 ? -1.0 : 1.0) * 0.0168077);
			}
			else { return sinusoidal(ipos); }
		}
		plnpos interruptedhomolosine(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd(), cmd;
			//the central meridian for each parts
			if (cyl.z > 0) {
				if (cyl.x < -40.0 * rdpdg) { cmd = spcpos(-100.0, 0.0, 0.0) * rdpdg; }
				else { cmd = spcpos(30.0, 0.0, 0.0) * rdpdg; }
			}
			else {
				if (cyl.x < -100.0 *rdpdg) { cmd = spcpos(-160.0, 0.0, 0.0) * rdpdg; }
				else if (cyl.x < -20.0 * rdpdg) { cmd = spcpos(-60.0, 0.0, 0.0) * rdpdg; }
				else if (cyl.x < 80.0 * rdpdg) { cmd = spcpos(20.0, 0.0, 0.0) * rdpdg; }
				else { cmd = spcpos(140.0, 0.0, 0.0) * rdpdg; }
			}
			return homolosine(mtrx3(cmd) * ipos) + plnpos(cmd.x / pi, 0.0);
		}
		plnpos aitoff(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd(), tpos = spcpos(sph.x / 2.0, sph.y, 1.0).XYsphpos();
			plnpos opos = azimuthalequidistant(tpos); opos.x *= 2.0;
			return opos;
		}
		plnpos hammer(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd(), tpos = spcpos(sph.x / 2.0, sph.y, 1.0).XYsphpos();
			plnpos opos = azimuthalequalarea(tpos) * sqrt(0.5); opos.x *= 2.0;
			return opos;
		}
		plnpos americanpolyconic(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			if (ipos.z == 0) { return plnpos(sph.x, 0) / pi; }
			else {
				double e = sph.x * ipos.z;
				return plnpos(sin(e) / tan(sph.y), sph.y + (1.0 - cos(e)) / tan(sph.y)) / pi;
			}
		}
		plnpos rectangularpolyconic(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd(); double a;
			if (dflat.x == 0) { a = sph.x / 2.0; }
			else { a = tan(sph.x * sin(dflat.x) / 2.0) / sin(dflat.x); }
			if (ipos.z == 0) { return plnpos(a * 2.0, 0) / fctr; }
			else {
				double e = atan(a * ipos.z) * 2.0;
				return plnpos(sin(e) / tan(sph.y), sph.y + (1.0 - cos(e)) / tan(sph.y)) / fctr;
			}
		}
		plnpos lagrange(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd();
			double w = dflat.y / pi, a1 = exp(log(2.0 / (1.0 - sin(dflat.x)) - 1.0) / w / 2.0);
			double f = tan(pi / w / 2.0) * 2.0;
			if (abs(cyl.z) == 1.0) { return plnpos(0.0, (cyl.z > 0 ? 2.0 : -2.0)) / f; }
			else {
				double a = exp(log(2.0 / (1.0 - cyl.z) - 1.0) / w / 2.0), v = a / a1;
				double c = (v + 1 / v) / 2.0 + cos(cyl.x / w);
				return plnpos(sin(cyl.x / w) * 2.0, v - 1 / v) / c / f;
			}
		}
		plnpos eisenlohr(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			double s1 = sin(sph.x / 2.0), c1 = cos(sph.x / 2.0);
			double t = sin(sph.y / 2.0) / (cos(sph.y / 2.0) + sqrt(cos(sph.y) * 2.0) * c1);
			double c = sqrt(2.0 / (1 + t * t));
			double v = sqrt((cos(sph.y / 2.0) + sqrt(cos(sph.y) / 2.0) * (c1 + s1)) / (cos(sph.y / 2.0) + sqrt(cos(sph.y) / 2.0) *(c1 - s1)));
			return plnpos(-2.0 * log(v) + c * (v - 1.0 / v), -2.0 * atan(t) + c * t * (v + 1.0 / v)) *0.935;
		}
		plnpos augustepicycloidal(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			double t2 = tan(sph.y / 2.0);
			double c1 = sqrt(1.0 - t2 * t2), c = 1.0 + c1 * cos(sph.x / 2.0);
			double x1 = sin(sph.x / 2.0) * c1 / c, y1 = tan(sph.y / 2.0) / c;
			return plnpos(x1 * (3.0 + x1 * x1 - y1 * y1 * 3.0), y1 * (3.0 + x1 * x1 * 3.0 - y1 * y1)) / 4.0;
		}
		plnpos eckert1(spcpos ipos) { spcpos sph = ipos.XYsphcrd(); return plnpos(sph.x * (1.0 - abs(sph.y) / pi), sph.y) / pi; }
		plnpos eckert2(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd(); double t = sqrt(1.0 - 0.75 * abs(cyl.z)), s = (cyl.z >= 0 ? 1.0 : -1.0);
			return plnpos(cyl.x * t / pi, s * (1.0 - t));
		}
		plnpos eckert3(spcpos ipos) { spcpos sph = ipos.XYsphcrd(); return plnpos((0.5 + sqrt(0.25 - sph.y * sph.y / pi / pi)) * sph.x, sph.y) / pi; }
		plnpos eckert4(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd(); double y = (pi + 4.0) * cyl.z / 2.0;
			double t = inversefnc(y, cyl.z / 2.5, [](double s) { return sin(s * 2.0) / 2.0 + sin(s) * 2.0 + s; }, [](double s) { return cos(s * 2.0) + cos(s) * 2.0 + 1.0; });
			return plnpos(cyl.x * (1.0 + cos(t)) / pi, sin(t)) / 2.0;
		}
		plnpos eckert5(spcpos ipos) { spcpos sph = ipos.XYsphcrd(); return plnpos((1.0 + cos(sph.y)) / 2.0 * sph.x, sph.y) / pi; }
		plnpos eckert6(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd(); double y = (pi / 2.0 + 1.0) * cyl.z;
			double t = inversefnc(y, cyl.z / 1.5, [](double s) { return sin(s) + s; }, [](double s) { return cos(s) + 1.0; });
			return plnpos(cyl.x * (1.0 + cos(t)) / 2.0, t) / pi;
		}
		plnpos kavrayskiy7(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			return plnpos(sph.x * 3.0 / pi * sqrt(pi * pi / 3.0 - sph.y * sph.y), sph.y * 2.0) / pi / sqrt(3.0);
		}
		plnpos craigretroazimuthal(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd();
			if (abs(cyl.x) < exceptionalmargin) { return plnpos(0, cyl.z - tan(dflat.x) * cyl.y) / pi; }
			else { return plnpos(cyl.x, cyl.x / sin(cyl.x) * (cyl.z * cos(cyl.x) - tan(dflat.x) * cyl.y)) / pi; }
		}
		plnpos hammerretroazimuthal(spcpos ipos) {
			spcpos cyl = ipos.XYcylcrd();
			double csr = sin(dflat.x) * ipos.z + cos(dflat.x) * ipos.x;
			double cstsnr = cos(dflat.x) * ipos.y / cyl.y;
			double sntsnr = -sin(dflat.x) * cyl.y + cos(dflat.x) * ipos.z * ipos.x / cyl.y;
			double snr = sqrt(cstsnr * cstsnr + sntsnr * sntsnr);
			return plnpos(cstsnr, sntsnr) / snr * atan2(snr, csr) / pi;
		}
		plnpos littrow(spcpos ipos) {
			double rr = ipos.x * ipos.x + ipos.y * ipos.y;
			return plnpos(ipos.y, ipos.x * ipos.z) / rr / pi;
		}
		plnpos loximuthal(spcpos ipos) {
			spcpos sph = ipos.XYsphcrd();
			if (abs(sph.y - dflat.x) < exceptionalmargin) { return plnpos(sph.x * cos(sph.y), sph.y) / pi; }
			else { return plnpos(sph.x * (sph.y - dflat.x) / log(tan(pi / 4.0 + sph.y / 2.0) / tan(pi / 4.0 + dflat.x / 2.0)), sph.y) / pi; }
		}
		plnpos peircequincuncial(spcpos ipos) {//integral of an original definition, so very slow.
			int nm = 256;//integral step
			plnpos z = plnpos(ipos.y, ipos.z) / (ipos.x + 1.0);
			if (z.sqlength() == 0) { return plnpos(0, 0); }
			if (ipos.x < 0) { z = plnpos(1.0, 0.0) / z; }//for good convergence
			double r = sqrt(z.sqlength()), stp = 1.0 / (double)nm;
			plnpos rt, rt2, ss, z2, zr, dz = z * (stp / r), zz = dz / 2.0;
			int ch = 0;
			for (double i = stp / 2.0; i < r; i += stp) {//integral
				z2 = zz * zz; z2 = (plnpos(1.0, 0.0) - z2 * z2); zr = z2.complexsqrt();
				ss += dz / zr;
				if (ch == 0) { if (z2.sqlength() < 0.6) { ch = 1; stp /= 2.0; dz /= 2.0; } }
				else if (ch == 1) { if (z2.sqlength() < 0.25) { ch = 2; stp /= 2.0; dz /= 2.0; } }
				else if (ch == 2) { if (z2.sqlength() < 0.09) { ch = 3; stp /= 2.0; dz /= 2.0; } }
				zz += dz;
			}
			if (ipos.x < 0) {
				if (ipos.y * ipos.z > 0) { rt = plnpos(0.0, 1.0); }
				else { rt = plnpos(1.0, 0.0); }
				if (ipos.y > 0) { rt *= -1.0; }
				rt2 = plnpos(0.0, 1.0) * rt;
				ss = (ss * rt + plnpos(1.31103, 1.31103)) * rt; ss = plnpos(0.0, 1.0) * ss;
			}
			return ss / 1.31103;
		}
		plnpos guyou(spcpos ipos) {
			plnpos opos = peircequincuncial(mtrx3::R045 * ipos);
			return plnpos(opos.x - opos.y, opos.x + opos.y);
		}
		plnpos test(spcpos ipos) {
			plnpos opos = stereographic(ipos);
			return opos.complexsqrt() * 0.5;
		}
};