/*=============================================================================
    acme_magfield.cpp	- get the magnetic field at a given point,
	                    also the BTransverse for a given velocity

	-------------------
    begin				: 05-22-04
    copyright			: (C) 2004 by Andrew O'Neill, HiRes Grad Student
    email				: oneill@phys.columbia.edu
=============================================================================*/
 /* This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program 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 General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */

#include "acme_magfield.h"

using namespace std;

Vector3D MagField (const Vector3D& Position)
{
   // GEOPACK uses geocentric polar coords, and float <-> real
   float R = Position.R()/ (GeopackEarthRadius + 0.0); // R must be in units of R_e
//   float Theta = Position.Lat(); // CHANGED 2005-03-21 from Theta() and Phi()
// 	float Phi = Position.Long();  
   float Theta = Position.Theta(); 
 	float Phi = Position.Phi();  

   // Create variables to store result of geopack in.
   // These are Outward, Southward and Eastward components 
   // (i.e. Theta, Phi are not angles!)  in nTesla
   float BR = 0.0, BTheta = 0.0, BPhi = 0.0;

   // Get the magnetic field from GEOPACK in nT
#ifdef INTEL
   igrf_geo_(R, Theta, Phi, BR, BTheta, BPhi);
#else
   igrf_geo__(R, Theta, Phi, BR, BTheta, BPhi);
#endif
//	dout(LINFO) << "B(R, Theta, Phi) = (" << BR << ", " << BTheta << ", " 
//	           << BPhi << ")\n";
//	dout(LINFO) << "|B| = " << sqrt(BR*BR + BTheta*BTheta + BPhi*BPhi) << "\n";

   // convert this to our geocentric coordinate system
   // Note that BTheta is in the southwards direction... 
   // i.e. -BTheta = BNorth
 
   Vector3D B; 
   B.setENU(BPhi, -BTheta, BR, Position.Lat(), Position.Long());

	B = B*1E-9;   // convert to Tesla (GEOPACK uses nT)
//	dout(LINFO) << "B (x, y, z) = " << B << endl;
//	dout(LINFO) << "|B| = " << B.R() << endl;
	return B;
}

double TransverseMagField(const Vector3D& Position, const Vector3D& Velocity)
{
	Vector3D B = MagField(Position);
	double BTransverse =  B.Transverse(Velocity);
//	dout(LINFO) << "Position " << Position << " B_t: " << BTransverse << endl;
	return BTransverse;
}

void InitMagField (void)
{
	time_t ltime = time(NULL);
   struct tm *Now = localtime(&ltime);
	int IYear = Now->tm_year + 1900;
   int IDay = Now->tm_yday + 1;
	dout(INFO) << "Initialising mganetic field for Y" << IYear << " d" << IDay
	           << " " << Now->tm_hour << ":" << Now->tm_min << ":" << Now->tm_sec
	           << "\n";
   recalc_(IYear , IDay, Now->tm_hour, Now->tm_min, Now->tm_sec);
}