//Current version
///ID:          $Id: Coord.cs,v 1.1.1.1 2006/11/14 02:56:46 gupchup Exp $
///Revision:    $Revision: 1.1.1.1 $
///Date:        $Date: 2006/11/14 02:56:46 $

using System;
using System.Drawing;
using System.Drawing.Drawing2D;

namespace Geom
{
	/// <summary>
	/// Summary description for Coord.
	/// </summary>
	/// <summary>
	/// Coord private class carries the transformation of ra/dec to nu/mu and xy. 
	/// </summary>
	public class Coord 
	{
				
		public static string Revision
		{
			get 
			{
				return "$Revision: 1.1.1.1 $";
			}
		}
			
		public double ra,dec;			// point (ra, dec) in J2000 degrees
		public double mu,nu;			// point (ra, dec) in J2000 degrees
		public double row, col;			// the row and column number in the frame

		private float x,y;				// pixel offsets from upper left corner 
		// corresponding to (ra,dec), scaled 
		public double scale;			// image scale 
		private double a,b,c,d,e,f,node,incl;	// affine tranformation of the stripe to ra,dec
		public string info;				// information about the field for debugging
		public Matrix m;
		public static double D2R = Math.PI/180.0;		// degrees to radians

		
		public Coord()
		{}

		public Coord(double va,double vb,double vc,double vd,double ve,double vf,
			double vnode,double vincl,double vscale, string vinfo) 
		{
			a = va;
			b = vb;
			c = vc;
			d = vd;
			e = ve;
			f = vf;
			node = vnode;
			incl = vincl;
			scale= vscale;
			info = vinfo;
			m	 = null;
		}

		/// <summary>
		/// The copy method copies the values of the coord parameter into the current coordinate object
		/// </summary>
		/// <param name=coord>			
		public void copy( Coord coord ) 
		{
			this.a			= coord.a;
			this.b			= coord.b;
			this.c			= coord.c;
			this.d			= coord.d;
			this.e			= coord.e;
			this.f			= coord.f;
			this.node		= coord.node;
			this.incl		= coord.incl;
			this.scale		= coord.scale;
			this.info		= coord.info;
			this.m			= coord.m;
		}

		/// <summary>
		/// EqtoScreen goes from equatorial coordinates to xy
		/// Sets ra, dec and consequently xOffset and yOffset of 
		/// Coord according to the astrometric transformation.
		/// </summary>
		/// <param name="ra">right ascension degrees in J2000</param>
		/// <param name="dec">declination degrees in J2000</param>
		public PointF EqToFrame(double pRa, double pDec) 
		{   
			// save the angles
			ra	= pRa;
			dec	= pDec;

			// convert angles to radians
			double inclR = incl*D2R;
			double raR	 = (pRa-node) *D2R;
			double decR	 = pDec*D2R;

			// first go to (mu,nu)
			double gx	=  Math.Cos(raR) * Math.Cos(decR);
			double gy	=  Math.Sin(raR) * Math.Cos(decR);
			double gz	=  Math.Sin(decR); 	

			// reverse rotation by incl around x
			double qx	=  gx;
			double qy	=  gy*Math.Cos(inclR) + gz * Math.Sin(inclR);
			double qz	= -gy*Math.Sin(inclR) + gz * Math.Cos(inclR);

			// compute mu, nu
			mu	=  Math.Atan2(qy,qx)/D2R + node;
			nu	=  Math.Asin(qz)/D2R;

			// adjust for wraparound of the difference
			double dmu	=  mu - a;
			if (dmu<-180)  dmu += 360;
			if (dmu> 180)  dmu -=360; 
	
			// set up the determinant
			double det	 = b * f - e * c;
			double col   = ( (nu-d)*b - dmu*e )/det;
			double row   = ( dmu*f - (nu-d)*c )/det;
				
			// compute pixel coordinates
			x	= (float)(col/scale);
			y	= (float)(row/scale);

			PointF p = new PointF(x,y);
			return p; 
		} 

		/// <summary>
		/// FrameToEq goes from Frame x,y coordinates and sets ra and dec
		/// </summary>
		/// <param name="x">pixels across (scaled for zoom level)</param>
		/// <param name="y">pixels down (scaled for zoom level)</param>
		public PointEq FrameToEq(float x, float y)  
		{
			// x,y need to be in original pixel scale
			double col = x;
			double row = y;
    
			// Great circle coordinates in degrees
			mu = (a +  b*row  +  c*col);
			nu = (d +  e*row  +  f*col);

			// convert to radians, then to ra,dec
			double muR = (mu - node) * D2R;
			double nuR = nu   * D2R;
			double xg  = Math.Cos(muR) * Math.Cos(nuR);
			double yg  = Math.Sin(muR) * Math.Cos(nuR);
			double zg  = Math.Sin(nuR);

			// compute unit vector
			double inR = incl * D2R;
			double qx   = xg;
			double qy   = yg * Math.Cos(inR) - zg*Math.Sin(inR);
			double qz   = yg * Math.Sin(inR) + zg*Math.Cos(inR);

			// compute ra, dec
			double pRa  = Math.Atan2(qy,qx) / D2R  +  node;
			double pDec = Math.Asin(qz) / D2R;
			if (pRa>360) pRa -= 360;
			if (pRa<0)   pRa += 360;

			ra	= pRa;
			dec	= pDec;

			return new PointEq(pRa,pDec);
		}

		/// <summary>
		/// rotateDegreesToNorth returns the number of degrees 
		/// the image must be rotated to make North be "up"
		/// </summary>
		/// <param name="ra"> center point of the frame</param>
		/// <param name="dec">center point of the frame</param>
		public double rotateDegreesToNorth(double ra, double dec)  
		{
			// compute point due North, and center
			PointF n = this.EqToFrame(ra,dec+0.1);
			PointF c = this.EqToFrame(ra,dec);			
			double angle = -90 - Math.Atan2((double)(n.Y-c.Y),(double)(n.X-c.X)) /D2R;
			return angle;
		}
	} // end of Coord Class
}
/* Revision History

        $Log: Coord.cs,v $
        Revision 1.1.1.1  2006/11/14 02:56:46  gupchup
        no message

        Revision 1.1.1.1  2006/04/25 16:07:23  gupchup
        Code base that merges original imgcutout code and WCS frame format imgcutout code

        Revision 1.3  2003/04/08 17:12:00  nieto
        Trying to include projection.cs class in project

        Revision 1.2  2003/03/25 23:18:33  nieto
        testing draw fields with their own coord
        
		
*/
