// GeometricAlgebra Java package file
   // Current version
   // Go to ... http://sinai.mech.fukui-u.ac.jp/gcj/software/KamiWaAi/index.html
   // Copyright (C) 2003, Eckhard M.S. Hitzer
   //
   // This library is free software; you can redistribute it and/or
   // modify it under the terms of the GNU Lesser General Public
   // License as published by the Free Software Foundation; either
   // version 2.1 of the License, or any later version.
  //
  // This library is distributed in the hope that it will be useful,
  // but WITHOUT ANY WARRANTY; without even the implied warranty of
  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  // Lesser General Public License for more details.
  //
  // You should have received a copy of the GNU Lesser General Public
  // License along with this library; if not, write to the Free Software
  // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  // Go to ... http://www.gnu.org/copyleft/lesser.html
  // 
  // hitzer@mech.fukui-u.ac.jp
  // 
  // Department of Mechanical Engineering, Fukui University
  // 3-9-1 Bunkyo, 910-8507 Fukui, Japan
  // 
  //
  //
  // It the personal wish of E.M.S. Hitzer,
  // that you never use this software for infringing human dignity:
  //
  // "God created man in his own image,
  // in the image of God he created him;
  // male and female he created them.
  // [The Bible, Genesis chapter 1 verse 27]
  //
  //
  //
  // 
  //File Sphere.java
  // 2003-5-1
  // radius(), Center(), center() methods tested on 2003-1-29
  package net.sourceforge.kamiwaai.geometricalgebra;
  import java.awt.Color;
  import java.awt.Graphics;
  import java.util.Vector;
  public class Sphere extends GeometricObject {
  	private MultiVector[] basepoints = new MultiVector[4];
  	// The default color of circles will be magenta.
  	private Color color = Color.magenta;
  	public Sphere(MultiVector L4) {
  		L = L4;
  	}
  	public Sphere(PointC point1, PointC point2, PointC point3, PointC point4) {
  		L = point1.getMultiVector().OutProd(point2.getMultiVector()).OutProd(
  				point3.getMultiVector().OutProd(point4.getMultiVector()));
  	}
  	public Sphere(PointC centre, double radius){
  		L = I.mult(centre.getMultiVector().sub(n.multSc(0.5*radius*radius))).multSc(2.0*radius*radius*radius);
  	}
  	public double radius() {
  		double r2 = (L.ScProd(L)) * (1.0 / (n.OutProd(L).ScProd(n.OutProd(L))));
  		double r = Math.sqrt(r2);
  		return r;
  	}
  	public MultiVector Center() {
  		// Gives the full conformal MultiVector for the center of the sphere
  		MultiVector low = n.OutProd(L).mult(I).multSc(-1.0);
  		double llow = low.ScProd(low);
  		MultiVector Low = low.multSc(1.0 / llow); // inverse of low
  		MultiVector Result = L.mult(I).mult(Low);
  		double r = radius();
  		MultiVector Cent = Result.add(n.multSc(0.5 * r * r));
  		return Cent;
  	}
  	public double[] center() {
  		// Gives the 3D Euclidean vector of the center of the circle
  		double[] cvector = new double[3];
  		Sphere LS = new Sphere(L);
  		cvector[0] = LS.Center().get3DEuclidianPoint()[0];
  		cvector[1] = LS.Center().get3DEuclidianPoint()[1];
  		cvector[2] = LS.Center().get3DEuclidianPoint()[2];
  		return cvector;
  	}
  	public MultiVector generator(int inc, MultiVector plane) {
  		// Gives the Motor (MuliVector that rotates a given
  		// angle increment 2pi/inc about the sphere's center)
  		int UP = inc;
  		double angle = 2 * Math.PI / UP; 
  		MultiVector cv, Iplane, Tc, Tcr, R, D;
  		cv = Center().OutProd(N).mult(N); // Euclidean 3D center vector.
  		// Translation operator
  		Tc = One.add(n.mult(cv).multSc(0.5));
  		Tcr = Tc.reverse();
  		Iplane = plane.multSc(1.0 / (plane.magnitude()));
  		// Rotation rotor R:
  		R = One.multSc(Math.cos(0.5 * angle)).add(
  				Iplane.multSc(Math.sin(0.5 * angle)));
  		// Operator for rotation about the circle center by angle increment.
  		D = Tc.mult(R).mult(Tcr);
 		return D;
 	}
 	public MultiVector PonS() {
 		// Gives back one point on the sphere, which is e.g. used as the Pole
 		// for generating a long. and lat. net of circles on the sphere,
 		// and for picking the sphere itself.
 		MultiVector rv, rad, Trad, Tradr, pons;
 		double r;
 		// Maybe later rv will be an input parameter for PonS
 		rv = e1.add(e2).add(e3.multSc(-1.0));
 		r = radius();
 		rad = rv.multSc(r / (rv.magnitude()));
 		// Translation operator from the center to a point on the cirlce itself
 		Trad = One.add(n.mult(rad).multSc(0.5));
 		Tradr = Trad.reverse();
 		pons = Trad.mult(Center()).mult(Tradr);
 		return pons;
 	}
 	public MultiVector PonSsouth() {
 		// Gives back the point on the sphere opposite to the PonS() point Pole.
 		// This makes the picking of a sphere more comfortable
 		MultiVector rv, rad, Trad, Tradr, pons;
 		double r;
 		// Maybe later rv will be an input parameter for PonS
 		rv = e1.add(e2).add(e3.multSc(-1.0));
 		r = radius();
 		rad = rv.multSc(r / (rv.magnitude()));
 		// Translation operator from the center to a point on the cirlce itself
 		Trad = One.add(n.mult(rad).multSc(0.5));
 		Tradr = Trad.reverse();
 		pons = Tradr.mult(Center()).mult(Trad);
 		// The translation of the Center is opposite to the one in PonS()
 		return pons;
 	}
 	public Vector snet(int inc) {
 		// Gives back list of lat. and long. circles on the sphere
 		Vector scnet = new Vector();
 		MultiVector pole, center, A1, A2, A3, L3, D1, D1r, D2, D2r;
 		MultiVector Iequp, Dlg, Dlgr;
 		Circle C, C0;
 		//Sphere LS = new Sphere(L);
 		pole = PonS().OutProd(N).mult(N);
 		center = Center().OutProd(N).mult(N);
 		D1 = generator(2 * inc, pole.sub(center).OutProd(e1));
 		D1r = D1.reverse();
 		D2 = generator(2 * inc, pole.sub(center).OutProd(e2));
 		D2r = D2.reverse();
 		// Producing the latitudinal meridians:
 		A1 = PonS();
 		A2 = PonS();
 		A3 = PonS();
 		for (int k = 1; k < inc; k++) {
 			A1 = D1.mult(A1).mult(D1r);
 			A2 = D2.mult(A2).mult(D2r);
 			A3 = D1r.mult(A3).mult(D1);
 			L3 = A1.OutProd(A2).OutProd(A3);
 			C = new Circle(L3);
 			C.setColor(color);
 			//if(k==1) Iequp = C.Cplane();
 			scnet.add(C);
 		}
 		// This defines a 3D Euclidean bivector parallel to the
 		// plane of the equator on the sphere relative to the Pole.
 		Iequp = ((Circle) scnet.get(0)).Cplane();
 		// Producing the longitudinal circles:
 		A1 = PonS();
 		A2 = D1.mult(A1).mult(D1r);
 		A3 = D1r.mult(A1).mult(D1);
 		L3 = A1.OutProd(A2).OutProd(A3);
 		C0 = new Circle(L3);
 		C0.setColor(color);
 		scnet.add(C0);
 		Dlg = generator(2 * inc, Iequp);
 		Dlgr = Dlg.reverse();
 		for (int k = 1; k < inc; k++) {
 			C = new Circle((Dlg.Powerof(k)).mult(C0.getMultiVector()).mult(
 					Dlgr.Powerof(k)));
 			C.setColor(color);
 			scnet.add(C);
 		}
 		return scnet;
 		// Description of snet() method:
 		// Define a set of inc horizontal (latitude) planes hplanes
 		// Define a set of inc vertical (longitude) planes all containing rv
 		// Which is longitudinal, which lateral? Name accordingly
 		// Intersect these planes with the sphere LS to get the
 		// net of circles (a vector of circles) which can be drawn with
 		// the drawC() and drawCZ() methods of the circle objects.
 		// Test every method with a simple example in GeoPro
 	}
 	public void setColor(Color c) {
 		color = c;
 	}
 	public void draw(Graphics g) {
 		Sphere LS = new Sphere(L);
 		Vector SNET;
 		SNET = LS.snet(6);
 		for (int k = 0; k <= SNET.size() - 1; k++) {
 			((Circle) SNET.get(k)).setColor(color);
 			((Circle) SNET.get(k)).draw(g);
 		}
 	}
 	public void drawZ(Graphics g) {
 		Sphere LS = new Sphere(L);
 		Vector SNET;
 		SNET = LS.snet(6);
 		for (int k = 0; k <= SNET.size() - 1; k++) {
 			((Circle) SNET.get(k)).setColor(color);
 			((Circle) SNET.get(k)).drawZ(g);
 		}
 	}
 	public void setBasePoint(MultiVector basepoint, int p1to4) {
 		basepoints[p1to4] = basepoint;
 	}
 	// This method is specific for remaking a sphere when one of its
 	// basepoints is shifted!
 	public void remake(MultiVector movedBP, int bpno) {
 		basepoints[bpno] = movedBP;
 		MultiVector Lnew = basepoints[0].OutProd(basepoints[1]).OutProd(
 				basepoints[2].OutProd(basepoints[3]));
 		L = Lnew;
 	}
 	public void recenter(MultiVector NewCenter) {
 		MultiVector OldCenter, cDiff, Tc, Tcr, Lnew;
 		Sphere LS = new Sphere(L);
 		OldCenter = LS.Center();
 		cDiff = NewCenter.get3DMVector().sub(OldCenter.get3DMVector());
 		Tc = One.add(n.mult(cDiff).multSc(0.5));
 		Tcr = Tc.reverse();
 		Lnew = Tc.mult(L).mult(Tcr);
 		L = Lnew;
 	}
 	private int MeetNo(Line line) {
 		int meetno;
 		double meet2 = meet(line).ScProd(meet(line));
 		if (meet2 > 0.0)
 			meetno = 2;
 		else if (meet2 == 0.0)
 			meetno = 1;
 		else
 			meetno = 0;
 		return meetno;
 	}
 	public Vector MeetLineMVs(Line line) {
 		Vector mline = new Vector();
 		//Sphere S = new Sphere(L);
 		int mno = MeetNo(line);
 		//mline.add(mno); // This would add the number of intersection points
 		MultiVector Meet, Meetp, u, v, a, b;
 		double u2, gamma;
 		//double[] a3D = new double[3];
 		//double[] b3D = new double[3];
 		Meet = meet(line);
 		Meetp = Meet.sub(Meet.mult(N).getGrade(0).mult(N)); // Chopping off the
 															// N part
 		v = Meetp.mult(nbar).multSc(-1.0).getGrade(1).multSc(-2.0);
 		u = Meetp.mult(n).multSc(-1.0).getGrade(1);
 		u2 = u.ScProd(u);
 		gamma = 2.0 * Meet.getnhnbPart().getScPart().RealPart();
 		if (mno == 1) {
 			MultiVector fact, factinv, A;
 			fact = Meet.mult(n).multSc(-1.0).getGrade(1);
 			factinv = fact.multSc(1.0 / (fact.magnitude() * fact.magnitude()));
 			// This defines the full conformal vector aa:
 			A = factinv.mult(Meet);
 			// The spatial part of a seems to have the
 			// right form for the mno = 2 case "stem" part
 			//a3D[0] = -a.getnhnbPart().getBvPart()[0].ImaginaryPart();
 			//a3D[1] = -a.getnhnbPart().getBvPart()[1].ImaginaryPart();
 			//a3D[2] = -a.getnhnbPart().getBvPart()[2].ImaginaryPart();
 			//mline.add(a3D);
 			mline.add(A);
 		} else {
 			if (mno == 2 && gamma == 0.0) {
 				MultiVector Bv, udB, x, stem, A, B;
 				double ra = radius();
 				Bv = Meetp.mult(N).getGrade(4).mult(N);
 				udB = u.Lcontract(Bv);
 				stem = udB.multSc(1.0 / u2);
 				x = u.multSc(Math.sqrt(ra * ra / u2 - udB.ScProd(udB)
 						/ (u2 * u2 * u2)));
 				a = stem.add(x);
 				b = stem.sub(x);
 				A = a.add(n.multSc(0.5 * a.ScProd(a))).add(nbar);
 				B = b.add(n.multSc(0.5 * b.ScProd(b))).add(nbar);
 				//a3D[0] = -a.getnhnbPart().getBvPart()[0].ImaginaryPart();
 				//a3D[1] = -a.getnhnbPart().getBvPart()[1].ImaginaryPart();
 				//a3D[2] = -a.getnhnbPart().getBvPart()[2].ImaginaryPart();
 				//b3D[0] = -b.getnhnbPart().getBvPart()[0].ImaginaryPart();
 				//b3D[1] = -b.getnhnbPart().getBvPart()[1].ImaginaryPart();
 				//b3D[2] = -b.getnhnbPart().getBvPart()[2].ImaginaryPart();
 				mline.add(A);
 				mline.add(B);
 			} else if (mno == 2 && gamma != 0.0) {
 				MultiVector A, B;
 				double sigma, v2, rho, a2, b2;
 				sigma = 0.5 * gamma * gamma - u.ScProd(v);
 				v2 = v.ScProd(v);
 				rho = Math.sqrt(sigma * sigma - v2 * u2);
 				// In the following two lines, the sigma part produces
 				// the "stem" part
 				a2 = (sigma + rho) / u2;
 				b2 = (sigma - rho) / u2;
 				a = (u.multSc(a2).add(v)).multSc(1.0 / gamma);
 				b = (u.multSc(b2).add(v)).multSc(1.0 / gamma);
 				A = a.add(n.multSc(0.5 * a2)).add(nbar);
 				B = b.add(n.multSc(0.5 * b2)).add(nbar);
 				//a3D[0] = -a.getnhnbPart().getBvPart()[0].ImaginaryPart();
 				//a3D[1] = -a.getnhnbPart().getBvPart()[1].ImaginaryPart();
 				//a3D[2] = -a.getnhnbPart().getBvPart()[2].ImaginaryPart();
 				//b3D[0] = -b.getnhnbPart().getBvPart()[0].ImaginaryPart();
 				//b3D[1] = -b.getnhnbPart().getBvPart()[1].ImaginaryPart();
 				//b3D[2] = -b.getnhnbPart().getBvPart()[2].ImaginaryPart();
 				mline.add(A);
 				mline.add(B);
 			}
 		}
 		return mline;
 	}
 
 }