/*@@
   @file      prim2con_MHD.c
   @date      Sep 17, 2016
   @author    Philipp Moesta, Daniel Siegel
   @desc 
   prim2con for MHD 
   @enddesc 
 @@*/



#include <assert.h>
#include "con2prim.h"

extern struct c2p_steer c2p;

void prim2con_MHD_(const double * prim, double * con, const double g_con[4][4],
      const double g_cov[4][4]){
  
  
  // FORMAT
  // prim - rho, v_j, eps, B^k, Ye, T
  // con - D, S_j, tau, B^k, Ye_con, T

  // below is code for option to recompute Lorentz factor based on velocity, currently unused 

  // Raise indices
  //const double v1_con = g_con[1][1]*prim[v1_cov]+g_con[1][2]*prim[v2_cov]+g_con[1][3]*prim[v3_cov];
  //const double v2_con = g_con[2][1]*prim[v1_cov]+g_con[2][2]*prim[v2_cov]+g_con[2][3]*prim[v3_cov];
  //const double v3_con = g_con[3][1]*prim[v1_cov]+g_con[3][2]*prim[v2_cov]+g_con[3][3]*prim[v3_cov];

  // v^2 = v^i * v_i
  //const double v_squared = prim[v1_cov]*v1_con + prim[v2_cov]*v2_con + prim[v3_cov]*v3_con;
  
  // Lorentz factor
  //const double W = 1.0/sqrt(1.0-v_squared);

  const double W = prim[WLORENTZ];
  assert(W >= 1.0);
 
  // Anton et al. Equation (44)
  const double alpha_b0 = W*(prim[B1_con]*prim[v1_cov]+prim[B2_con]*prim[v2_cov]+prim[B3_con]*prim[v3_cov]);
  
  // Lower indices - covariant
  const double B1_cov = g_cov[1][1]*prim[B1_con]+g_cov[1][2]*prim[B2_con]+g_cov[1][3]*prim[B3_con];
  const double B2_cov = g_cov[2][1]*prim[B1_con]+g_cov[2][2]*prim[B2_con]+g_cov[2][3]*prim[B3_con];
  const double B3_cov = g_cov[3][1]*prim[B1_con]+g_cov[3][2]*prim[B2_con]+g_cov[3][3]*prim[B3_con];
  
  // Anton et al. Equations (44) & (45)
  const double b1_cov = B1_cov/W + alpha_b0*prim[v1_cov];
  const double b2_cov = B2_cov/W + alpha_b0*prim[v2_cov];
  const double b3_cov = B3_cov/W + alpha_b0*prim[v3_cov];
  
  // B^2 = B^i * B_i
  const double B_squared = prim[B1_con]*B1_cov + prim[B2_con]*B2_cov + prim[B3_con]*B3_cov;
  
  // Anton et al. Equation (46)
  const double b_squared = (B_squared + pow(alpha_b0,2))/pow(W,2);
 
  //double P = P_prim_(c2p.eoskey,prim);
  const double P_star = prim[PRESS] + b_squared/2.0;
  const double eps_star = prim[EPS] + b_squared/(2.0*prim[RHO]);
  double hstar;
  if(c2p.eoskey==1){ // Polytropic EOS, h = h(rho)
    hstar = 1.0 + b_squared/(2.0*prim[RHO]) + P_star/prim[RHO];
  } else { // h = h(rho, eps)
    hstar = 1.0 + eps_star + P_star/prim[RHO];
  }
  
  // Anton et al. Equation (37)
  con[D] = prim[RHO]*W;
  
  // Anton et al. Equation (38)
  // S_j = rho*h*W^2*v_j - a*b^0*b_j
  // rho*h = rho + rho*eps + press
  con[S1_cov] = prim[RHO]*hstar*prim[v1_cov]*pow(W,2) - alpha_b0*b1_cov;
  con[S2_cov] = prim[RHO]*hstar*prim[v2_cov]*pow(W,2) - alpha_b0*b2_cov;
  con[S3_cov] = prim[RHO]*hstar*prim[v3_cov]*pow(W,2) - alpha_b0*b3_cov;
  
  // Anton et al. Equation (39)
  // tau = rho*h*W^2 - p - (a*b^0)^2 - D
  con[TAU] = prim[RHO]*hstar*pow(W,2) - P_star - pow(alpha_b0,2) - con[D]; 
  // + con[D] to recover Noble & Gammie
  
  con[B1_con] = prim[B1_con];
  con[B2_con] = prim[B2_con];
  con[B3_con] = prim[B3_con];
  
  con[YE] = prim[YE]*con[D];
  con[TEMP] = prim[TEMP];

}