#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include <unistd.h>
#include <getopt.h>

#define PI M_PI 
#define R0 8.3 //The distance of the Sun to the Galactic Center

// Galactic cartesian coordinates (x,y,z)
// x parallel to (l,b)(90,0)
// y parallel to (l,b)(180,0)
// GC (0,0,0)
// x_sun=0 kpc
// y_sun=8.3 kpc

// Galactic cylindrical coordinates (r,phi,z)
// GC(0,0,0)
// phi starts from x in CCW (x parallel to l=90)

//////////////functions
void diskB(double x, double y, double z, double *Br, double *Bphi, double *Bz);
//disk magnetic field modeled by Han et al. (2018)
//Br: radial component / Positive: outward
//Bphi: azimuthal component / Positive: CCW
//Bz: vertical component / Positive: north

void haloB(double x, double y, double z, double *Br, double *Bphi, double *Bz);
//halo toroidal magnetic field modeled by Xu & Han (2024)

void GMF(double l, double b, double d, double *Bt, double *phi, double *Bsi);
//Galactic magnetic field at any point (l,b,d)
//Bt: total amplitude of the field
//phi: azimuthal angle of the field measured counterclockwise from the x-axis
//Bsi: B component along the line of sight to the Sun / Positive: toward us

//////////////main function
int main(int argc,char *argv[])
{

  double l,b,d;
  l=0;b=0;d=2;
  
  if(argc<7)
    {
      printf("Usage: ./GMFcal -l 10 -b 10 -d 2 \n");
      return 0;
    }
  int opt,index;
  const char *optstring="l:b:d:";//here ':' denotes to have a parameter
  struct option opts[]={
    {"help",         no_argument ,   NULL,'h'},
    {   0,           0,              0,    0 }};

  while((opt=getopt_long(argc,argv,optstring,opts,&index))!=-1)
    {
      switch(opt)    
	{  	  
	case 'l':
	  l=atof(optarg);
	  break;
	case 'b':
	  b=atof(optarg);
	  break;
	case 'd':
	  d=atof(optarg);
	  break;
	default:       
	  break;
	}
    }

  if(l<0 || l>360)
    {
      printf("ERROR: Please input 0 <= l < 360 !\n");
      return 0;
    }
  if(b<-90 || b>90)
    {
      printf("ERROR: Please input -90 <= b <= 90 !\n");
      return 0;
    }


  l=l*PI/180.0;//rad
  b=b*PI/180.0;//rad

  double x,y,z;
  double phixy;
  double rxy;//Galactocentric radius

  double dxy;//distance to sun in the Galactic plane

  double Br=0;//radial component positive for outward direction
  double Bphi=0;//azimuthal component positive for counter-clockwise
  double Bz=0;

  //-disk field
  double Brd=0;
  double Bphid=0;
  double Bzd=0;
  
  //-halo toroidal field
  double Brt=0;
  double Bphit=0;
  double Bzt=0;
  
  dxy=d*cos(b);
  x=dxy*sin(l);
  y=R0-dxy*cos(l);
  z=d*sin(b);
  rxy=sqrt(x*x+y*y);
  
  if(rxy<=20 && rxy>0)
    {
      diskB(x, y, z, &Brd, &Bphid, &Bzd);
      haloB(x, y, z, &Brt, &Bphit, &Bzt);
      
      Br=Brd+Brt;
      Bphi=Bphid+Bphit;
      Bz=Bzd+Bzt;
    }
  else
    {
      Br=0;
      Bphi=0;
      Bz=0;
    }


  printf("Disk magnetic field (microG): Brd=%.2f Bphid=%.2f Bzd=%0.2f \n", Brd,Bphid,Bzd);
  printf("Halo magnetic field (microG): Brh=%.2f Bphih=%.2f Bzh=%0.2f \n", Brt,Bphit,Bzt);
  printf("Total magnetic field (microG): Br=%.2f Bphi=%.2f Bz=%0.2f \n", Br,Bphi,Bz);
  
}

void diskB(double x, double y, double z, double *Br, double *Bphi, double *Bz)
{
  //parameters initialization
  double p=11.0 * PI/180.0; //pitch angle in rad
  double Hd=0.4; //vertical scalehight
  double rd=5.0; //radial scalelength
  double Rr[]={3.0, 4.1, 4.9, 6.1, 7.5, 8.16, 10.5, 15.0};//Radial Zones
  double Bs[]={4.5,-3.0, 6.3, -4.7, 3.3, -8.7, 4.5};//CCW: positive; CW:negative

  //variable parameters
  double Bt=0; //total disk B
  double rxy; //Galactocentric radius
  double phixy; //azimutal angle measured counterclockwise from +y-axis
  double r0; //the radius on the positive Y-axis from the Galactic Centre

  double tanp;
  tanp=tan(p);
  
  rxy=sqrt(x*x+y*y);

  if(rxy<3.0 || rxy>20.0) //out of region set Bt=0
    {
      *Bphi=0;
      *Br=0;
    }
  else
    {
      phixy=atan2(y,x)*180./PI-90.; //from +y-axis
      if(phixy <= 0) phixy=phixy+360;

      r0=rxy/exp(tanp*(phixy)/180*PI);
      if(r0>=Rr[0] && r0<Rr[1]) Bt=Bs[0]*exp(-rxy/rd); 
      if(r0>=Rr[1] && r0<Rr[2]) Bt=Bs[1]*exp(-rxy/rd); 
      if(r0>=Rr[2] && r0<Rr[3]) Bt=Bs[2]*exp(-rxy/rd); 
      if(r0>=Rr[3] && r0<Rr[4]) Bt=Bs[3]*exp(-rxy/rd); 
      if(r0>=Rr[4] && r0<Rr[5]) Bt=Bs[4]*exp(-rxy/rd); 
      if(r0>=Rr[5] && r0<Rr[6]) Bt=Bs[5]*exp(-rxy/rd); 
      if(r0>=Rr[6] && r0<Rr[7]) Bt=Bs[6]*exp(-rxy/rd);

      if(r0<=Rr[0])
	{ r0=rxy/exp(tanp* (phixy-360)/180*PI);
	  if(r0>=Rr[0] && r0<Rr[1]) Bt=Bs[0]*exp(-rxy/rd); 		    
	  if(r0>=Rr[1] && r0<Rr[2]) Bt=Bs[1]*exp(-rxy/rd); 
	  if(r0>=Rr[2] && r0<Rr[3]) Bt=Bs[2]*exp(-rxy/rd); 
	  if(r0>=Rr[3] && r0<Rr[4]) Bt=Bs[3]*exp(-rxy/rd); 
	  if(r0>=Rr[4] && r0<Rr[5]) Bt=Bs[4]*exp(-rxy/rd); 
	  if(r0>=Rr[5] && r0<Rr[6]) Bt=Bs[5]*exp(-rxy/rd); 
	  if(r0>=Rr[6] && r0<Rr[7]) Bt=Bs[6]*exp(-rxy/rd); 
	}

      if(r0>=Rr[6])
	{ r0=rxy/exp(tanp* (phixy+360)/180*PI);
	  if(r0>=Rr[0] && r0<Rr[1]) Bt=Bs[0]*exp(-rxy/rd);	  
	  if(r0>=Rr[1] && r0<Rr[2]) Bt=Bs[1]*exp(-rxy/rd);
	  if(r0>=Rr[2] && r0<Rr[3]) Bt=Bs[2]*exp(-rxy/rd); 
	  if(r0>=Rr[3] && r0<Rr[4]) Bt=Bs[3]*exp(-rxy/rd); 
	  if(r0>=Rr[4] && r0<Rr[5]) Bt=Bs[4]*exp(-rxy/rd); 
	  if(r0>=Rr[5] && r0<Rr[6]) Bt=Bs[5]*exp(-rxy/rd); 
	  if(r0>=Rr[6] && r0<Rr[7]) Bt=Bs[6]*exp(-rxy/rd); 
	}
      
      *Bphi=Bt*cos(p);
      *Br=Bt*sin(p);
    }
  *Bphi=*Bphi *exp(-fabs(z)/Hd);
  *Br=*Br *exp(-fabs(z)/Hd);
  *Bz=0;
}

void haloB(double x, double y, double z, double *Br, double *Bphi, double *Bz)
{
  //parameters initialization
  double z0=3.00; //vertical scale height in kpc
  double B0=0.73; //maximum field strength at Rt and Rh in muG
  double Rt=5.31; //scale radius in kpc
  double Rh=7.97; //in kpc
  
  //variable parameters
  double rxy;
  rxy=sqrt(x*x+y*y);

  double haloZ,haloR;
  double halo_field;

  double sign(double x) //sign(z) function for z>0 positive, counter-clockwise
  { 
    if(x>0) 
      return 1; 
    else if(x<0)
      return -1;
    else 
      return 0; 
  }
  
  haloZ=(fabs(z)/z0)*exp(-(fabs(z)-z0)/z0);
  haloR=exp(-(rxy-Rh)/Rt*(rxy-Rh)/Rt);
  
  halo_field=sign(z)*B0*haloZ*haloR;
      
  *Br=0;
  *Bphi=halo_field;
  *Bz=0;
}

void GMF(double l, double b, double d, double *Bt, double *phi, double *Bsi)
//input l,b in degree
//input d in kpc
//"+" toward us;  "-" opposite
{
  l=l*PI/180.0;//rad
  b=b*PI/180.0;//rad

  double x,y,z;
  double phixy;
  double rxy;//Galactocentric radius

  double dxy;//distance to sun in the Galactic plane

  double Br=0;//radial component positive for outward direction
  double Bphi=0;//azimuthal component positive for counter-clockwise
  double Bz=0;

  //-disk field
  double Brd=0;
  double Bphid=0;
  double Bzd=0;
  
  //-halo toroidal field
  double Brt=0;
  double Bphit=0;
  double Bzt=0;
  
  dxy=d*cos(b);
  x=dxy*sin(l);
  y=R0-dxy*cos(l);
  z=d*sin(b);
  rxy=sqrt(x*x+y*y);
  
  if(rxy<=20 && rxy>0)
    {
      diskB(x, y, z, &Brd, &Bphid, &Bzd);
      haloB(x, y, z, &Brt, &Bphit, &Bzt);
      
      Br=Brd+Brt;
      Bphi=Bphid+Bphit;
      Bz=Bzd+Bzt;
    }
  else
    {
      Br=0;
      Bphi=0;
      Bz=0;
    }
  
  //B component in the sight of line
  double angle_dr;//the angle between dxy and rxy
  
  phixy=atan2(y,x);//-pi--pi
  if(phixy<0)
    phixy=phixy+2*PI;//pi--2pi  
  angle_dr=(phixy+PI/2)-l;//triangle
  *Bsi=                   //"+" toward us
    (-Br*cos(angle_dr)  //r component
     +Bphi*sin(angle_dr))//phi component
    *cos(b)
    -Bz*sin(b);//z component

  
  //total B
  *Bt=sqrt(Br*Br+Bphi*Bphi+Bz*Bz);
  //azimuthal angle of B
  double phib;
  phib=atan2(Bphi,Br);//-pi--pi
  if(phib<0)
    phib=phib+2*PI;//pi--2pi
  *phi=(phib+phixy)*180/PI; //in degree
  if(*phi>360)
    *phi=*phi-360;

  if((rxy*rxy+z*z)>(20*20) || rxy==0) //if out of 20 kpc from the GC set 0
    {
      *Bt=0;
      *phi=0;
      *Bsi=0;
    }
}