text:
#if !defined(THDM_H)
#define THDM_H

#include "SM.h"
#include <complex>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix.h>

#include "CouplingFunctionTypeI.h"
#include "CouplingFunctionTypeII.h"
#include "CouplingFunctionTypeLS.h"
#include "CouplingFunctionTypeFL.h"

using namespace std;

/**
* \mainpage 2HDMC Class documentation
* 2HDMC is a general-purpose calculator for the two-Higgs doublet model.
* It allows parametrization of the Higgs potential in many different ways,
* convenient specification of generic Yukawa sectors, the evaluation of decay
* widths, theoretical constraints and much more.
*/

/**
* @brief Specifies a general two-Higgs doublet model
*
* This class implements a description of a general two-Higgs doublet model
* in terms of the potential parameters and Yukawa couplings necessary for
* a full specification of the model. Several methods are available to set
* the model parameters to various special cases.
*
* From a THDM object, all the couplings between physical states can be
* accessed, including Higgs-fermion, Higgs-Vector, Higgs-Higgs, 3H and 4H
* couplings.
*/
//static int thdmc_set_param(int key, int npara, double *para, double *res);

int thdmc_set_param(int key, double smpara[], double para[], double res[], int slha);

typedef ComplexType (*CFunc_NLO)(double, double, double, double, double, double, double, double, double, double, SM);
ComplexType ZEROCOUPLING(double MHL2, double MHH2, double MHA2, double MHp2, double M2, double beta, double alpha, double m12, double m22, double m32, SM sm);
// For simplicity, even though SVV, SSV and SSS coupling is independent of the Yukawa type, I still use them to indicate the formula.
#define UPID 10
#define DOID 20
#define LEID 30
#define NUID 40
enum NLOSTATUS
{
    NLOGOOD = 1,
    NOTIMPLEMENT = 0
};
int GetNLOFuncSFF(int S1,int F1,int F2,int type,CFunc_NLO &func);
int GetNLOFuncSFFC(int S1,int F1,int F2,int type,CFunc_NLO &func_PL, CFunc_NLO &func_PR);
int GetNLOFuncSVV(int S1,int V1,int V2,int type,CFunc_NLO &func_gmunu, CFunc_NLO &func_kmuqnu);
int GetNLOFuncSSV(int S1,int S2,int V1,int type,CFunc_NLO &func);
int GetNLOFuncSSS(int S1,int S2,int S3,int type,CFunc_NLO &func);

class THDM {

 public:

  /**
  * @brief Default constructor.
  *
  * Empty default constructor which does nothing except initialization of this
  * THDM object. The model is not specified in any way by using this constructor.
  */
  THDM();
  // ~THDM();

  /**
  * @brief Sets the underlying SM
  *
  * This method can be used to specify the SM parameters which this 2HDM
  * is using.
  *
  * @param sm_in SM object specifying the parameters to use in this 2HDM
  */
  void set_SM(SM sm_in);

  /**
  * @brief Returns the underlying SM
  *
  * @returns The SM object specifying the parameters used in this 2HDM
  */
  SM get_SM();

  /**
  * @brief Specifies 2HDM using generic potential
  *
  * This method lets the user specify the 2HDM using a generic basis.
  * For details on the basis choices available for the Higgs potential, we
  * refer to the complete manual.
  *
  * @param lambda1 Value of \f$ \lambda_1 \f$
  * @param lambda2 Value of \f$ \lambda_2 \f$
  * @param lambda3 Value of \f$ \lambda_3 \f$
  * @param lambda4 Value of \f$ \lambda_4 \f$
  * @param lambda5 Value of \f$ \lambda_5 \f$
  * @param lambda6 Value of \f$ \lambda_6 \f$
  * @param lambda7 Value of \f$ \lambda_7 \f$
  * @param m12_2 Soft \f$ Z_2 \f$-breaking parameter
  * @param tan_beta Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_gen(double lambda1, double lambda2, double lambda3,
                     double lambda4, double lambda5, double lambda6,
                     double lambda7, double m12_2, double tan_beta);

  /**
  * @brief Specifies 2HDM in the Higgs basis
  *
  * This method lets the user specify the 2HDM in the Higgs basis \f$ (v_2\equiv 0) \f$.
  * For details on the basis choices available for the Higgs potential, we
  * refer to the complete manual.
  *
  * @param Lambda1 Value of \f$ \Lambda_1 \f$
  * @param Lambda2 Value of \f$ \Lambda_2 \f$
  * @param Lambda3 Value of \f$ \Lambda_3 \f$
  * @param Lambda4 Value of \f$ \Lambda_4 \f$
  * @param Lambda5 Value of \f$ \Lambda_5 \f$
  * @param Lambda6 Value of \f$ \Lambda_6 \f$
  * @param Lambda7 Value of \f$ \Lambda_7 \f$
  * @param m_Hp    Mass of the charged Higgs
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_higgs(double Lambda1, double Lambda2, double Lambda3,
                       double Lambda4, double Lambda5, double Lambda6,
                       double Lambda7, double m_Hp);

  /**
  * @brief Specifies 2HDM in the hybrid basis of 1507.04281
  *
  * @param mh Value of light CP-even Higgs mass \f$ m_h \f$
  * @param mH Value of heavy CP-even Higgs mass \f$ m_H \f$
  * @param cba Mixing parameter \f$ \cos(\beta-\alpha) \f$. Valid range: -1 < cba <= 1.
  * @param Z_4 Value of \f$ Z_4 \f$
  * @param Z_5 Value of \f$ Z_5 \f$
  * @param Z_7 Value of \f$ Z_7 \f$
  * @param tanb Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_hybrid(double mh, double mH, double cba, double Z4,double Z5, double Z7, double tanb);

  /**
  * @brief Specifies 2HDM in the hybrid basis of 1507.04281, but with sba instead of cba
  *
  * @param mh Value of light CP-even Higgs mass \f$ m_h \f$
  * @param mH Value of heavy CP-even Higgs mass \f$ m_H \f$
  * @param sba Mixing parameter \f$ \sin(\beta-\alpha) \f$. Valid range: -1 < cba <= 1.
  * @param Z_4 Value of \f$ Z_4 \f$
  * @param Z_5 Value of \f$ Z_5 \f$
  * @param Z_7 Value of \f$ Z_7 \f$
  * @param tanb Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_hybrid_sba(double mh, double mH, double sba, double Z4,double Z5, double Z7, double tanb);

  /**
  * @brief Specifies 2HDM using potential from the Higgs Hunter's Guide
  *
  * This method lets the user specify the 2HDM using the CP-conserving form
  * of the potential given by Eq. (4.8) in "The Higgs Hunter's Guide".
  * For details on the basis choices available for the Higgs potential, we
  * refer to the complete manual.
  *
  * @param lambda1 Value of \f$ \lambda_1 \f$
  * @param lambda2 Value of \f$ \lambda_2 \f$
  * @param lambda3 Value of \f$ \lambda_3 \f$
  * @param lambda4 Value of \f$ \lambda_4 \f$
  * @param lambda5 Value of \f$ \lambda_5 \f$
  * @param lambda6 Value of \f$ \lambda_6 \f$
  * @param tan_beta Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_HHG(double lambda1, double lambda2, double lambda3,
                     double lambda4, double lambda5, double lambda6,
                     double tan_beta);

  /**
  * @brief Specifies 2HDM in the physical basis
  *
  * This method lets the user specify the 2HDM using the basis of physical
  * Higgs masses. For details on the basis choices available for the Higgs
  * potential, we refer to the complete manual.
  *
  * @param m_h  Mass of lightest CP-even Higgs \f$ h \f$
  * @param m_H  Mass of heavier CP-even Higgs \f$ H \f$
  * @param m_A  Mass of CP-odd Higgs \f$ A \f$
  * @param m_Hp Mass of charged Higgs
  * @param sba  Mixing parameter \f$ \sin(\beta-\alpha) \f$. NB: Correct sign on \f$ \sin(\beta-\alpha) \f$ must be determined from the condition \f$ \cos(\beta-\alpha) \geq 0 \f$
  * @param lambda6  Value of \f$ \lambda_6 \f$ in generic potential
  * @param lambda7  Value of \f$ \lambda_7 \f$ in generic potential
  * @param m12_2    Soft \f$ Z_2 \f$-breaking parameter
  * @param tan_beta Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_param_phys(double m_h,double m_H, double m_A, double m_Hp,
                      double sba, double lambda6, double lambda7,
                      double m12_2, double tan_beta);

  bool set_param_sm(double mh);

  /**
  * @brief Specifies the 2HDM of the tree-level MSSM
  *
  * This method lets the user specify a 2HDM with the properties of the tree-level
  * MSSM in terms of masses and coupling relations. The Yukawa sector is also
  * automatically selected to be of type II.
  *
  * @param m_A      Mass of CP-odd Higgs \f$ A \f$
  * @param tan_beta Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_MSSM(double m_A, double tan_beta);

  /**
  * @brief Specifies the 2HDM of the tree-level MSSM + mass corrections to h ("hMSSM")
  *
  * This method lets the user specify a 2HDM with the properties of the tree-level
  * MSSM, with additional mass corrections from the lambda_2 contribution to the (2,2)
  * element of the mass matrix. These can be identified as the leading (mt^4) MSSM
  * corrections. The Yukawa sector is automatically selected to be of type II.
  *
  * @param m_h      Mass of light CP-even Higgs \f$ h \f$
  * @param m_A      Mass of CP-odd Higgs \f$ A \f$
  * @param tan_beta Ratio of vevs, \f$ \tan\beta=v_2/v_1 \f$
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_hMSSM(double mh, double mA, double tanb);

  /**
  * @brief Specifies the 2HDM with one "inert Higgs" doublet
  *
  * This method lets the user specify a 2HDM with special properties as follows:
  * only one of the doublets acquires a vev, there is an exact \f$ Z_2 \f$-symmetry,
  * preventing interdoublet mixing, and the doublet without a vev has no Yukawa
  * couplings. This results in one doublet being "inert" (or dark). The lightest
  * Higgs originating from this doublet is then stable, thus a dark matter candidate.
  * The conventions used for specifying the parameters for this model are the same as
  * in arXiv:0810.3924.
  *
  * @param m_h  Mass of "SM-like" Higgs - NB: Changed meaning compared to "physical" basis
  * @param m_H  Mass of CP-even "inert" Higgs - NB: Changed meaning compared to "physical" basis
  * @param m_A  Mass of CP-odd "inert" Higgs
  * @param m_Hp Mass of charged "inert" Higgs
  * @param lambda2 Value of \f$ \lambda_2 \f$ parameter
  * @param lambda3 Value of \f$ \lambda_2 \f$ parameter
  *
  * @returns Boolean \a true if all parameters were set correctly, \a false otherwise
  */
  bool set_inert(double m_h,double m_H, double m_A, double m_Hp, double lambda2, double lambda3);

  /**
  * @brief Returns parameter set in the generic basis
  *
  * This method returns a consistent set of parameter values describing the
  * current model in the generic basis.
  *
  * @param lambda1 Returned value of \f$ \lambda_1 \f$
  * @param lambda2 Returned value of \f$ \lambda_2 \f$
  * @param lambda3 Returned value of \f$ \lambda_3 \f$
  * @param lambda4 Returned value of \f$ \lambda_4 \f$
  * @param lambda5 Returned value of \f$ \lambda_5 \f$
  * @param lambda6 Returned value of \f$ \lambda_6 \f$
  * @param lambda7 Returned value of \f$ \lambda_7 \f$
  * @param m12_2   Returned value of \f$ m_{12}^2 \f$
  * @param tan_beta Returned value of \f$ \tan\beta \f$
  */
  void get_param_gen(double &lambda1, double &lambda2, double &lambda3,
                     double &lambda4, double &lambda5, double &lambda6,
                     double &lambda7, double &m12_2,   double &tan_beta);

  /**
  * @brief Returns parameter set in the Higgs basis
  *
  * This method returns a consistent set of parameter values describing the
  * current model in the Higgs basis.
  *
  * @param Lambda1 Returned value of \f$ \Lambda_1 \f$
  * @param Lambda2 Returned value of \f$ \Lambda_2 \f$
  * @param Lambda3 Returned value of \f$ \Lambda_3 \f$
  * @param Lambda4 Returned value of \f$ \Lambda_4 \f$
  * @param Lambda5 Returned value of \f$ \Lambda_5 \f$
  * @param Lambda6 Returned value of \f$ \Lambda_6 \f$
  * @param Lambda7 Returned value of \f$ \Lambda_7 \f$
  * @param m_Hp    Returned value of the charged Higgs mass
  */
  void get_param_higgs(double &Lambda1, double &Lambda2, double &Lambda3,
                       double &Lambda4, double &Lambda5, double &Lambda6,
                       double &Lambda7, double &m_Hp);

  void get_param_hybrid(double &m_h, double &m_H, double &sba,
                        double &Z4, double &Z5, double &Z7, double &tan_beta);

  /**
  * @brief Returns parameter set in Higgs Hunter's Guide basis
  *
  * This method returns a consistent set of parameter values describing the
  * current model in the basis used in the Higgs Hunter's Guide. NB: The use
  * of this method assumes \f$ \lambda_6=\lambda_7=0 \f$ in the generic basis.
  *
  * @param lambda1 Returned value of \f$ \lambda_1 \f$
  * @param lambda2 Returned value of \f$ \lambda_2 \f$
  * @param lambda3 Returned value of \f$ \lambda_3 \f$
  * @param lambda4 Returned value of \f$ \lambda_4 \f$
  * @param lambda5 Returned value of \f$ \lambda_5 \f$
  * @param lambda6 Returned value of \f$ \lambda_6 \f$
  * @param tan_beta Returned value of \f$ \tan\beta \f$
  */
  void get_param_HHG(double &lambda1, double &lambda2, double &lambda3,
                     double &lambda4, double &lambda5, double &lambda6,
                     double &tan_beta);

  /**
  * @brief Returns parameter set in physical basis
  *
  * This method returns a consistent set of parameter values describing the
  * current model in the "physical" basis.
  *
  * @param m_h      Returned value of lightest CP-even Higgs mass
  * @param m_H      Returned value of heaviest CP-even Higgs mass
  * @param m_A      Returned value of CP-odd Higgs mass
  * @param m_Hp     Returned value of charged Higgs mass
  * @param sba      Returned value of \f$ \sin(\beta-\alpha) \f$
  * @param lambda6  Returned value of \f$ \lambda_6 \f$
  * @param lambda7  Returned value of \f$ \lambda_7 \f$
  * @param m12_2    Returned value of \f$ m_{12}^2 \f$
  * @param tan_beta Returned value of \f$ \tan\beta \f$
  */
  void get_param_phys(double &m_h,double &m_H, double &m_A, double &m_Hp,
                      double &sba, double &lambda6, double &lambda7,
                      double &m12_2, double &tan_beta);

  /**
  * @brief Changes basis for the Higgs doublets
  *
  * This method performs a change of basis, changing the value of
  * \f$ \tan\beta \f$. This results in a recalculation of all the other
  * potential parameters, but the physical couplings remain unchanged.
  * Yukawa couplings are not modified.
  *
  * @param tan_beta Value of \f$ \tan\beta \f$ to use in the new basis
  */
  void recalc_tan_beta(double tan_beta);

  /**
  * @brief Returns the mass of physical Higgs bosons
  *
  * This method returns the mass of a physical Higgs Boson. The numbering
  * convention corresponds to (1,2,3,4 = h,H,A,H+) where always \f$ m_h < m_H \f$.
  *
  * @param h Index of Higgs boson
  *
  * @returns The mass of Higgs boson \a h
  */
  double get_hmass(int h);

  /**
  * @brief Returns invariant \f$ \sin(\beta-\alpha) \f$
  *
  * @returns Value of \f$ \sin(\beta-\alpha) \f$ in model
  */
  double get_sba();

  /**
  * @brief Returns invariant \f$ \cos(\beta-\alpha) \f$
  *
  * @returns Value of \f$ \cos(\beta-\alpha) \f$ in model
  */
  double get_cba();

  /**
  * @brief Higgs coupling \f$ q_{ki} \f$ factors
  *
  * This method returns the invariants \f$ q_{ki} \f$ (hep-ph/0602242) which are used
  * for the triple and quartic Higgs couplings.
  *
  * @param k Higgs index (1--4)
  * @param i Coupling index (1--2)
  *
  * @returns The value of the coefficient \f$ q_{ki} \f$
  */
  complex <double> get_qki(int k, int i);

  /**
  * @brief Initializes complete Yukawa sector to a specific type
  *
  * This method is used to specify a type of Yukawa sector for the 2HDM.
  * The types (1-4) implemented follow the convention of .
  *
  * @param type Type of Yukawa sector (1--4)
  */
  void set_yukawas_type(int type);

  /**
  * @brief Initializes and sets diagonal Yukawa couplings for down-type quarks
  *
  * This method initializes and sets the Yukawa couplings for the down-type quarks.
  * The diagonal matrix \f$ \kappa^D \f$ is automatically specified to comply with
  * the quark pole masses, whereas the diagonal elements of \f$ \rho^D \f$ can be
  * specified as input.
  *
  * @param rhod Diagonal element of \f$ \rho^D \f$ for the \f$ d \f$ quark
  * @param rhos Diagonal element of \f$ \rho^D \f$ for the \f$ s \f$ quark
  * @param rhob Diagonal element of \f$ \rho^D \f$ for the \f$ b \f$ quark
  *
  * @see set_yukawas_up, set_yukawas_lepton
  */
  void set_yukawas_down(double rhod, double rhos, double rhob);

  /**
  * @brief Initializes and sets diagonal Yukawa couplings for up-type quarks
  *
  * This method initializes and sets the Yukawa couplings for the up-type quarks.
  * The diagonal matrix \f$ \kappa^U \f$ is automatically specified to comply with
  * the quark pole masses, whereas the diagonal elements of \f$ \rho^U \f$ can be
  * specified as input.
  *
  * @param rhou Diagonal element of \f$ \rho^U \f$ for the \f$ u \f$ quark
  * @param rhoc Diagonal element of \f$ \rho^U \f$ for the \f$ c \f$ quark
  * @param rhot Diagonal element of \f$ \rho^U \f$ for the \f$ t \f$ quark
  *
  * @see set_yukawas_down, set_yukawas_lepton
  */
  void set_yukawas_up(double rhou, double rhoc, double rhot);

  /**
  * @brief Initializes and sets diagonal Yukawa couplings for the charged leptons
  *
  * This method initializes and sets the Yukawa couplings for the leptons.
  * The diagonal matrix \f$ \kappa^L \f$ is automatically specified to comply with
  * the lepton pole masses, whereas the diagonal element of \f$ \rho^L \f$ can be
  * specified as input.
  *
  * @param rhoe   Diagonal element of \f$ \rho^L \f$ for \f$ e \f$
  * @param rhomu  Diagonal element of \f$ \rho^L \f$ for \f$ \mu \f$
  * @param rhotau Diagonal element of \f$ \rho^L \f$ for \f$ \tau \f$
  *
  * @see set_yukawas_down, set_yukawas_up
  */
  void set_yukawas_lepton(double rhoe, double rhomu, double rhotau);

  /**
  * @brief Initializes and sets Yukawa couplings for down-type quarks
  *
  * This method initializes and sets the Yukawa couplings for the down-type quarks.
  * The diagonal matrix \f$ \kappa^D \f$ is automatically specified to comply with
  * the quark pole masses, whereas the symmetric matrix \f$ \rho^D \f$ is given by
  * the six elements specified as input.
  *
  * @param rho11 Yukawa coupling \f$ \rho^D_{dd} \f$
  * @param rho22 Yukawa coupling \f$ \rho^D_{ss} \f$
  * @param rho33 Yukawa coupling \f$ \rho^D_{bb} \f$
  * @param rho12 Yukawa coupling \f$ \rho^D_{ds} \f$
  * @param rho13 Yukawa coupling \f$ \rho^D_{db} \f$
  * @param rho23 Yukawa coupling \f$ \rho^D_{sb} \f$
  */
  void set_yukawas_down(double rho11,double rho22,double rho33,double rho12,double rho13,double rho23);

  /**
  * @brief Initializes and sets Yukawa couplings for up-type quarks
  *
  * This method initializes and sets the Yukawa couplings for the up-type quarks.
  * The diagonal matrix \f$ \kappa^U \f$ is automatically specified to comply with
  * the quark pole masses, whereas the symmetric matrix \f$ \rho^U \f$ is given by
  * the six elements specified as input.
  *
  * @param rho11 Yukawa coupling \f$ \rho^U_{uu} \f$
  * @param rho22 Yukawa coupling \f$ \rho^U_{cc} \f$
  * @param rho33 Yukawa coupling \f$ \rho^U_{tt} \f$
  * @param rho12 Yukawa coupling \f$ \rho^U_{uc} \f$
  * @param rho13 Yukawa coupling \f$ \rho^U_{ut} \f$
  * @param rho23 Yukawa coupling \f$ \rho^U_{ct} \f$
  */
  void set_yukawas_up(double rho11,double rho22,double rho33,double rho12,double rho13,double rho23);

  /**
  * @brief Initializes and sets Yukawa couplings for leptons
  *
  * This method initializes and sets the Yukawa couplings for the leptons.
  * The diagonal matrix \f$ \kappa^L \f$ is automatically specified to comply with
  * the lepton masses, whereas the symmetric matrix \f$ \rho^L \f$ is given by
  * the six elements specified as input.
  *
  * @param rho11 Yukawa coupling \f$ \rho^L_{ee} \f$
  * @param rho22 Yukawa coupling \f$ \rho^L_{\mu\mu} \f$
  * @param rho33 Yukawa coupling \f$ \rho^L_{\tau\tau} \f$
  * @param rho12 Yukawa coupling \f$ \rho^L_{e\mu} \f$
  * @param rho13 Yukawa coupling \f$ \rho^L_{e\tau} \f$
  * @param rho23 Yukawa coupling \f$ \rho^L_{\mu\tau} \f$
  */
  void set_yukawas_lepton(double rho11,double rho22,double rho33,double rho12,double rho13,double rho23);

  /**
  * @brief Initializes and sets Yukawa couplings for the inert 2HDM
  *
  * This method initializes and sets the Yukawa couplings for the "inert"
  * model. That means there are no Yukawa couplings of any type, except what
  * is necessary for the mass terms, i.e. all \f$ \rho = 0 \f$.
  */
  void set_yukawas_inert();

  /**
  * @brief Returns diagonal elements of Yukawa matrix for down-type quarks
  *
  * @param rhod Returned value of \f$ \rho^D_{dd} \f$
  * @param rhos Returned value of \f$ \rho^D_{ss} \f$
  * @param rhob Returned value of \f$ \rho^D_{bb} \f$
  */
  void get_yukawas_down(double &rhod, double &rhos, double &rhob);

  /**
  * @brief Returns diagonal elements of Yukawa matrix for up-type quarks
  *
  * @param rhou Returned value of \f$ \rho^U_{uu} \f$
  * @param rhoc Returned value of \f$ \rho^U_{cc} \f$
  * @param rhot Returned value of \f$ \rho^U_{tt} \f$
  */
  void get_yukawas_up(double &rhou, double &rhoc, double &rhot);

  /**
  * @brief Returns diagonal elements of Yukawa matrix for leptons
  *
  * @param rhoe Returned value of \f$ \rho^L_{ee} \f$
  * @param rhomu Returned value of \f$ \rho^L_{\mu\mu} \f$
  * @param rhotau Returned value of \f$ \rho^L_{\tau\tau} \f$
  */
  void get_yukawas_lepton(double &rhoe, double &rhomu, double &rhotau);

  void get_kappa_down(double &kd, double &ks, double &kb);
  void get_kappa_up(double &ku, double &kc, double &kt);
  void get_kappa_lepton(double &ke, double &kmu, double &ktau);

  void get_kappa_down(double mu, double &kd, double &ks, double &kb);
  void get_kappa_up(double mu, double &ku, double &kc, double &kt);
  void get_kappa_lepton(double mu, double &ke, double &kmu, double &ktau);
  void get_rho_down(double mu, double &rd, double &rs, double &rb);
  void get_rho_up(double mu, double &ru, double &rc, double &rt);
  void get_rho_lepton(double mu, double &re, double &rmu, double &rtau);

  /**
  * @brief Returns Yukawa matrix for down-type quarks
  *
  * @param rho_D_out Returned Yukawa matrix \f$ \rho^D \f$
  */
  void get_yukawas_down(gsl_matrix *rho_D_out);

  /**
  * @brief Returns Yukawa matrix for up-type quarks
  *
  * @param rho_U_out Returned Yukawa matrix \f$ \rho^U \f$
  */
  void get_yukawas_up(gsl_matrix *rho_U_out);

  /**
  * @brief Returns Yukawa matrix for leptons
  *
  * @param rho_L_out Returned Yukawa matrix \f$ \rho^L \f$
  */
  void get_yukawas_lepton(gsl_matrix *rho_L_out);

  /**
  * @brief Couplings of Higgses to down-type fermions
  *
  * Calculates the coupling \f$ hdd \f$ between one physical Higgs state, specified
  * by \a h, and two down-type quarks \a f1 and \a f2.
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param f1 First fermion (1,2,3 = d,s,b)
  * @param f2 Second fermion (1,2,3 = d,s,b)
  * @param cs Returned (complex) value for scalar coupling
  * @param cp Returned (complex) value for pseudoscalar coupling
  */
  void get_coupling_hdd(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);
  void get_coupling_hdd_NLO(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);

  /**
  * @brief Couplings of Higgses to up-type fermions
  *
  * Calculates the coupling \f$ huu \f$ between one physical Higgs state, specified
  * by \a h, and two up-type quarks \a f1 and \a f2.
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param f1 First fermion (1,2,3 = u,c,t)
  * @param f2 Second fermion (1,2,3 = u,c,t)
  * @param cs Returned (complex) value for scalar coupling
  * @param cp Returned (complex) value for pseudoscalar coupling
  */
  void get_coupling_huu(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);
  void get_coupling_huu_NLO(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);

  /**
  * @brief Couplings of Higgses to mixed type fermions
  *
  * Calculates the coupling \f$ hdu \f$ between one physical Higgs state, specified
  * by \a h (only relevant one is charged Higgs), and two quarks \a f1 and \a f2.
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param d Down-type fermion (1,2,3 = d,s,b)
  * @param u Up-type fermion (1,2,3 = u,c,t)
  * @param cs Returned (complex) value for scalar coupling
  * @param cp Returned (complex) value for pseudoscalar coupling
  */
  void get_coupling_hdu(int h,int d,int u,complex <double> &cs, complex <double> &cp);
  void get_coupling_hdu_NLO(int h,int d,int u,complex <double> &cs, complex <double> &cp);

  /**
  * @brief Couplings of Higgses to charged leptons
  *
  * Calculates the coupling \f$ hll \f$ between one physical Higgs state, specified
  * by \a h, and two leptons \a f1 and \a f2.
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param f1 First fermion (1,2,3 = \f$ e,\mu,\tau \f$)
  * @param f2 Second fermion (1,2,3 = \f$ e,\mu,\tau \f$)
  * @param cs Returned (complex) value for scalar coupling
  * @param cp Returned (complex) value for pseudoscalar coupling
  */
  void get_coupling_hll(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);
  void get_coupling_hll_NLO(int h,int f1,int f2,complex <double> &cs, complex <double> &cp);
  /**
  * @brief Couplings of Higgses to mixed leptons
  *
  * Calculates the coupling \f$ hl\nu_l \f$ between one physical Higgs state, specified
  * by \a h (only relevant one is charged Higgs), and two leptons \a f1 and \a f2.
  *
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param l  Charged lepton (1,2,3 = \f$ e,\mu,\tau \f$)
  * @param n  Neutrino (1,2,3 = \f$ \nu_e,\nu_\mu,\nu_\tau \f$)
  * @param cs Returned (complex) value for scalar coupling
  * @param cp Returned (complex) value for pseudoscalar coupling
  */
  void get_coupling_hln(int h,int l,int n,complex <double> &cs, complex <double> &cp);
  void get_coupling_hln_NLO(int h,int l,int n,complex <double> &cs, complex <double> &cp);

  /**
  * @brief Couplings of Higgses to pairs of vector bosons
  *
  * Calculates the coupling \f$ hV_1 V_2 \f$ between one physical Higgs state, specified
  * by \a h , and two vector bosons \a v1 and \a v2. NB: Neutral Higgses have no coupling
  * to photons implemented, but the loop mediated decays \f$ h\to \gamma \gamma \f$ can
  * nevertheless be calculated using the DecayTable class. Conventions are according to hep-ph/0602242.
  *
  * @param v1 Index of first vector boson (1,2,3 = \f$ \gamma,Z,W^+ \f$)
  * @param v2 Index of second vector boson (1,2,3 = \f$ \gamma,Z,W^+ \f$)
  * @param h  Index of Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param c  Returned (complex) value for coupling
  */
  void get_coupling_vvh(int v1,int v2,int h,complex <double> &c);
  void get_coupling_vvh_NLO(int v1,int v2,int h,complex <double> &c_gmunu,complex <double> &c_kmuqnu);

  /**
  * @brief Couplings of vector bosons to pairs of Higgses
  *
  * Calculates the coupling \f$ Vh_1 h_2 \f$ between two physical Higgs states, specified
  * by \a h1 and \a h2, and one vector boson \a v. Conventions are according to hep-ph/0602242.
  *
  * @param h1 Index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h2 Index of second Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param v  Index of vector boson (1,2,3 = \f$ \gamma,Z,W^+ \f$)
  * @param c  Returned (complex) value for coupling
  */
  void get_coupling_vhh(int v,int h1,int h2,complex <double> &c);
  void get_coupling_vhh_NLO(int v,int h1,int h2,complex <double> &c);

  /**
  * @brief Triple Higgs couplings
  *
  * Calculates the coupling \f$ h_1 h_2 h_3 \f$ between three physical Higgs states.
  * Conventions are according to hep-ph/0602242
  *
  * @param h1 Index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h2 Index of second Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h3 Index of third Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param c  Returned (complex) value for coupling
  */
  void get_coupling_hhh(int h1,int h2,int h3,complex <double> &c);
  void get_coupling_hhh_NLO(int h1,int h2,int h3,complex <double> &c);

  /**
  * @brief Couplings of two vector bosons and two Higgses
  *
  * Calculates the coupling \f$ V_1 V_2 h_1 h_2 \f$ between two vector bosons and
  * two physical Higgs states. Conventions are according to hep-ph/0602242
  *
  * @param v1 Index of first vector boson (1,2,3 = \f$ \gamma,Z,W^+ \f$)
  * @param v2 Index of second vector boson (1,2,3 = \f$ \gamma,Z,W^+ \f$)
  * @param h1 Index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h2 Index of second Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param c  Returned (complex) value for coupling
  */
  void get_coupling_vvhh(int v1,int v2,int h1,int h2,complex <double> &c);

  /**
  * @brief Quartic Higgs couplings
  *
  * Calculates the coupling \f$ h_1 h_2 h_3 h_4 \f$ between four physical Higgs states.
  * Conventions are according to hep-ph/0602242
  *
  * @param h1 Index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h2 Index of second Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h3 Index of third Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param h4 Index of fourth Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param c  Returned (complex) value for coupling
  */
  void get_coupling_hhhh(int h1,int h2,int h3,int h4,complex <double> &c);

  /**
  * @brief Calculates tree-level unitarity constraints
  *
  * This method calculates the eigenvalues of the S-matrix for Higgs-Higgs scattering
  * as defined in hep-ph/0508020.
  *
  * @returns Value of largest eigenvalue
  */
  double calc_unitarity();

  /**
  * @brief Checks tree-level unitarity constraints
  *
  * This method checks whether the parameters of the Higgs potential results
  * in an S-matrix for Higgs-Higgs scattering that fullfills tree-level unitarity,
  * as discussed in hep-ph/0508020.
  *
  * @param unitarity_limit optional value to be used for unitarity limit (default = \f$ 16\pi \f$)
  *
  * @returns Boolean \a true if the unitarity constraint is satisfied, \a false otherwise
  */
  bool check_unitarity(double unitarity_limit = 16*M_PI);

  /**
  * @brief Calculates largest quartic Higgs coupling
  *
  * This method calculates all the quartic Higgs boson couplings and returns the largest one as well as the indices of the four Higgs bosons
  *
  * @param gmax  Returned (complex) value for coupling
  * @param imax  Returned index of first Higgs boson (1,2,3,4 = h,H,A,H+)
  * @param jmax  Returned index of second boson (1,2,3,4 = h,H,A,H+)
  * @param kmax  Returned index of third boson (1,2,3,4 = h,H,A,H+)
  * @param lmax  Returned index of fourth boson (1,2,3,4 = h,H,A,H+)
  */
  void calc_perturbativity(complex <double> &gmax,int &imax,int &jmax,int &kmax,int &lmax);

  /**
  * @brief Checks perturbativity
  *
  * This method checks whether the couplings of the Higgs and Gauge bosons fullfills perturbativity
  *
  * @param perturbativity_limit optional value to be used for perturbativity limit (default = \f$ 4\pi \f$)
  *
  * @returns Boolean \a true if the perturbativity constraint is satisfied, \a false otherwise
  */
  bool check_perturbativity(double perturbativity_limit = 4*M_PI);

  /**
  * @brief Checks the stability of the Higgs potential
  *
  * This method checks whether the Higgs potential is stable, i.e. if no
  * directions exist in field space for which \f$ V\to -\infty \f$ for large
  * values of the fields. This is done through a combination of analytical and numerical
  * algorithms, depending on the generality of the specified potential.
  *
  * @returns Boolean \a true if the potential is deemed stable, \a false otherwise
  */
  bool check_stability();

  /**
  * @brief Prints the potential parameters in the generic basis to stdout
  */
  void print_param_gen();

  /**
  * @brief Prints the potential parameters in the Higgs basis to stdout
  */
  void print_param_higgs();

  /**
  * @brief Prints the potential parameters in the hybrid basis to stdout
  */
  void print_param_hybrid();

  /**
  * @brief Prints the potential parameters in the physical basis to stdout
  */
  void print_param_phys();

  /**
  * @brief Prints the potential parameters in the Higgs Hunter's Guide basis to stdout
  */
  void print_param_HHG();

  /**
  * @brief Prints the Yukawa matrices to stdout
  */
  void print_yukawas();

  void print_hdecay();

  /**
  * @brief Reads the 2HDM parameters from a LesHouches compliant file
  *
  * This method reads the parameters necessary to specify the 2HDM (and the %SM)
  * from a file complying with the LesHouches standard. There are several options
  * for how to prepare this file (see manual for details,
  * or the code distribution for examples).
  *
  * @param file The name of the file to read
  *
  * @returns Boolean \a true if the file could be opened and read correctly AND
  *          a 2HDM was correctly specified in the file, \a false otherwise
  */
  bool read_LesHouches(const char* file);

  /**
  * @brief Writes output in a LesHouches compliant file format
  *
  * This method writes the potential parameters and other selected output to
  * a specified file in LesHouches-compliant file format which can then be
  * further processed by other codes.
  *
  * @param file The name of the file which is to be written
  * @param fulldecay If \a true, decay modes of the Higgs bosons in LesHouches
  *                  format are written to the file
  * @param couplings If \a true, all couplings (> 200 values) for the Higgs bosons
  *                  are written to the file. Should be used when supplying input
  *                  for MadGraph/MadEvent 2HDMC model
  * @param qcd_on    Turns QCD corrections on or off, default is on. Should be
  *                  turned off when supplying input for MadGraph/MadEvent 2HDMC model
  */
  void write_LesHouches(const char* file, bool fulldecay, bool couplings, bool qcd_on=true, bool HBHS=true);

  void write_model(const char* file);

  double get_alpha();

  /**
	* @brief Small value
	*
	*  Minimum value used for widths, branching ratios etc. to determine when
	*  something should be considered zero.
	*/
  const static double EPS = 1E-15;

 private:
  double      lambda[8];
  double      beta;
  double      m22_2;
  double      sinba;
  bool        params_set;
  double      v2;
  gsl_matrix *kappa_D;
  gsl_matrix *kappa_U;
  gsl_matrix *kappa_L;
  gsl_matrix *rho_D;
  gsl_matrix *rho_U;
  gsl_matrix *rho_L;
  gsl_matrix *rho_N;
  int         yukawas_type;

  static bool first_run;
  const static char *version;

  SM sm;

  void init();
  double get_m12_2();
  void set_kappa();
  void set_kappa_D();
  void set_kappa_U();
  void set_kappa_L();

  void print_info();

// Added by Y.Wu
  bool Updated;
  bool HpmUpdated;
  double L_MHL2, L_MHH2, L_MHA2, L_MHP2;
  double L_MHL, L_MHH, L_MHA, L_MHP, L_sba, L_L6, L_L7, L_M122, L_tb;
  double L_beta, L_alp;
  double L_M2;
  gsl_matrix *rho_D_NLO;
  gsl_matrix *rho_U_NLO;
  gsl_matrix *RD_NLO;// = gsl_matrix_alloc(3,3);
  gsl_matrix *RU_NLO;// = gsl_matrix_alloc(3,3);

  // To save time, we can first check `Updated` to see if we need to rerun this function
  void store_param_phys();
  // void store_Hpmff_Couplings();

};

#endif


text:
/*
 *  projcl_run.c
 *  Magic Maps
 *
 *  Created by Evan Miller on 2/12/11.
 *  Copyright 2011 __MyCompanyName__. All rights reserved.
 *
 */

#include <projcl/projcl.h>
#include <projcl/projcl_warp.h>
#include "projcl_run.h"
#include "projcl_kernel.h"
#include "projcl_util.h"
#include "projcl_spheroid.h"

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#endif
#ifdef __linux__
#include <cblas.h>
typedef int __CLPK_integer;
typedef double __CLPK_doublereal;
extern int dgetrf_(__CLPK_integer *, __CLPK_integer *, __CLPK_doublereal *, __CLPK_integer *, __CLPK_integer *, __CLPK_integer *);
extern int dgetri_(__CLPK_integer *, __CLPK_doublereal *, __CLPK_integer *, __CLPK_integer *, __CLPK_doublereal *, __CLPK_integer *, __CLPK_integer *);
#endif

#define EPS7 1.e-7

#define SEC_TO_RAD 4.84813681109535993589914102357e-6

static cl_int pl_enqueue_kernel_albers_equal_area(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_american_polyconic(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_lambert_azimuthal_equal_area(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_lambert_conformal_conic(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_mercator(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_oblique_stereographic(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_robinson(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_transverse_mercator(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);
static cl_int pl_enqueue_kernel_winkel_tripel(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count);

struct pl_projection_info {
    PLProjection proj;
    char         name[80];
    cl_int (*func)(PLContext *pl_ctx,
                          cl_kernel kernel,
                          PLProjectionParams *params,
                          cl_mem xy_in,
                          cl_mem xy_out,
                          size_t count);
};

static struct pl_projection_info _pl_projection_info[] = {
    {
        .proj = PL_PROJECT_ALBERS_EQUAL_AREA,
        .name = "albers_equal_area",
        .func = &pl_enqueue_kernel_albers_equal_area,
    },
    {
        .proj = PL_PROJECT_AMERICAN_POLYCONIC,
        .name = "american_polyconic",
        .func = &pl_enqueue_kernel_american_polyconic,
    },
    {
        .proj = PL_PROJECT_LAMBERT_AZIMUTHAL_EQUAL_AREA,
        .name = "lambert_azimuthal_equal_area",
        .func = &pl_enqueue_kernel_lambert_azimuthal_equal_area,
    },
    {
        .proj = PL_PROJECT_LAMBERT_CONFORMAL_CONIC,
        .name = "lambert_conformal_conic",
        .func = &pl_enqueue_kernel_lambert_conformal_conic,
    },
    {
        .proj = PL_PROJECT_MERCATOR,
        .name = "mercator",
        .func = &pl_enqueue_kernel_mercator,
    },
    {
        .proj = PL_PROJECT_OBLIQUE_STEREOGRAPHIC,
        .name = "oblique_stereographic",
        .func = &pl_enqueue_kernel_oblique_stereographic,
    },
    {
        .proj = PL_PROJECT_ROBINSON,
        .name = "robinson",
        .func = &pl_enqueue_kernel_robinson,
    },
    {
        .proj = PL_PROJECT_TRANSVERSE_MERCATOR,
        .name = "transverse_mercator",
        .func = &pl_enqueue_kernel_transverse_mercator,
    },
    {
        .proj = PL_PROJECT_WINKEL_TRIPEL,
        .name = "winkel_tripel",
        .func = &pl_enqueue_kernel_winkel_tripel,
    }
};

struct pl_datum_info {
    double dx;
    double dy;
    double dz;
    double ex;
    double ey;
    double ez;
    double ppm;
};

/* Source: "WGS 84 Implementation Manual" */
static struct pl_datum_info pl_datum_params[] = {
    /*    Dx       Dy        Dz      Ex        Ey       Ez       m   */

    { /* PL_DATUM_WGS_84 */
          0.0,     0.0,    0.0,     0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_WGS_72 */
          0.0,     0.0,    4.5,     0.0,     0.0,    -0.554,  0.22    },
    { /* PL_DATUM_ED_50 */
        -87.0,   -98.0,  -121.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_ED_79 */
        -86.0,   -98.0,  -119.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_ED_87 */
        -82.5,   -91.7,  -117.7,    0.1338, -0.0625, -0.047,  0.045   },
    { /* PL_DATUM_AUSTRIA_NS */
        595.6,    87.3,   473.3,    4.7994,  0.0671,  5.7850, 2.555   },
    { /* PL_DATUM_BELGIUM_50 */
        -55.0,    49.0,  -158.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_BERNE_1873 */
        649.0,     9.0,   376.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_CH_1903 */
        660.1,    13.1,   369.2,    0.8048,  0.5777,  0.9522, 5.660   },
    { /* PL_DATUM_DANISH_GI_1934 */
        662.0,    18.0,   734.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_NOUV_TRIG_DE_FRANCE_GREENWICH */
       -168.0,   -60.0,   320.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_NOUV_TRIG_DE_FRANCE_PARIS */
       -168.0,   -60.0,   320.0,    0.0,     0.0,  8414.03,   0.0     },
    { /* PL_DATUM_POTSDAM */
        587.0,    16.0,   393.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_GGRS_87 */
        199.6,   -75.1,  -246.3,    0.0202,  0.0034,  0.0135, -0.015  },
    { /* PL_DATUM_HJORSEY_55 */
        -73.0,    46.0,   -86.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_IRELAND_65 */
        506.0,  -122.0,   611.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_ITALY_1940 */
       -133.0,   -50.0,    97.0,    0.0,     0.0, 44828.40,   0.0     },
    { /* PL_DATUM_NOUV_TRIG_DE_LUX */
       -262.0,    75.0,    25.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_NETHERLANDS_1921 */
        719.0,    47.0,   640.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_OSGB_36 */
        375.0,  -111.0,   431.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_PORTUGAL_DLX */
        504.1,  -220.9,   563.0,    0.0,     0.0,    -0.554,  0.220   },
    { /* PL_DATUM_PORTUGAL_1973 */
        227.0,    97.5,    35.4,    0.0,     0.0,    -0.554,  0.220   },
    { /* PL_DATUM_RNB_72 */
       -104.0,    80.0,   -75.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_RT_90 */
        424.3,   -80.5,   613.1,    4.3965, -1.9866,  5.1846, 0.0     },
    { /* PL_DATUM_NAD_27 */
         -8.0,   160.0,   176.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_NAD_83 */
          0.0,     0.0,     0.0,    0.0,     0.0,     0.0,    0.0     },
    { /* PL_DATUM_ETRS_89 */
          0.0,     0.0,     0.0,    0.0,     0.0,     0.0,    0.0     }
};

struct pl_matrix {
    double elements[4][4];
};

static struct pl_matrix pl_affine_transform_make(
        double Rx, double Ry, double Rz, double M,
        double Dx, double Dy, double Dz) {
    /* store in column-major order */
    struct pl_matrix matrix = {
        .elements = {
            {      M,  M * Rz,  -M * Ry, 0.f },
            {-M * Rz,       M,   M * Rx, 0.f },
            { M * Ry, -M * Rx,        M, 0.f },
            {     Dx,      Dy,       Dz, 1.f }
        }
    };
    return matrix;
}

static const char *_pl_proj_name(PLProjection proj) {
    int i=0;
    const char *name = NULL;
    for (i=0; i<sizeof(_pl_projection_info)/sizeof(_pl_projection_info[0]); i++) {
        if (_pl_projection_info[i].proj == proj) {
            name = _pl_projection_info[i].name;
            break;
        }
    }
    return name;
}

static double _pl_mlfn(double phi, double sphi, double cphi, double *en) {
	cphi *= sphi;
	sphi *= sphi;
	return(en[0] * phi - cphi * (en[1] + sphi*(en[2] + sphi*(en[3] + sphi*en[4]))));
}

static double _pl_qsfn(double sinphi, double e, double one_es) {
	double con = e * sinphi;
	return (one_es * (sinphi / (1.0 - con * con) - (.5 / e) * (log1p(-con) - log1p(con))));
}

static double _pl_msfn(double sinphi, double cosphi, double es) {
	return (cosphi / sqrt(1.0 - es * sinphi * sinphi));
}

static double _pl_tsfn(double phi, double sinphi, double e) {
	double con = e * sinphi;
	return (tan(.5 * (M_PI_2 - phi)) /
            pow((1.0 - con) / (1.0 + con), .5 * e));
}

static cl_int _pl_set_projection_kernel_args(PLContext *ctx, cl_kernel kernel,
        cl_mem xy_in, cl_mem xy_out, size_t count, PLSpheroidInfo *info, int *offset_ptr) {
	cl_int error = 0;
	cl_uint vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
	int offset = *offset_ptr;

	error |= pl_set_kernel_arg_mem(ctx, kernel, offset++, xy_in);
	error |= pl_set_kernel_arg_mem(ctx, kernel, offset++, xy_out);
	error |= pl_set_kernel_arg_uint(ctx, kernel, offset++, vec_count);

	if (!_pl_spheroid_is_spherical(info->tag)) {
        error |= pl_set_kernel_arg_float(ctx, kernel, offset++, info->ecc);
        error |= pl_set_kernel_arg_float(ctx, kernel, offset++, info->ecc2);
        error |= pl_set_kernel_arg_float(ctx, kernel, offset++, info->one_ecc2);
	}

	*offset_ptr = offset;

	return error;
}

cl_kernel pl_find_projection_kernel(PLContext *pl_ctx, PLProjection proj, int fwd, PLSpheroid ell) {
	char requested_name[128];
	if (fwd) {
		sprintf(requested_name, "pl_project_%s_%c", _pl_proj_name(proj), _pl_spheroid_is_spherical(ell) ? 's' : 'e');
	} else {
		sprintf(requested_name, "pl_unproject_%s_%c", _pl_proj_name(proj), _pl_spheroid_is_spherical(ell) ? 's' : 'e');
	}
	return pl_find_kernel(pl_ctx, requested_name);
}

static cl_int _pl_enqueue_kernel_1d(cl_command_queue queue, cl_kernel kernel, size_t dim_count) {
    cl_int error = CL_SUCCESS;
	error = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &dim_count, NULL, 0, NULL, NULL);
    return error;
}

cl_int pl_read_buffer(cl_command_queue queue, cl_mem xy_out_buf, float *xy_out, size_t out_count) {
    cl_int error;

    error = clEnqueueReadBuffer(queue, xy_out_buf, CL_TRUE, 0, out_count, xy_out, 0, NULL, NULL);
	if (error != CL_SUCCESS)
		return error;

	error = clFinish(queue);
	if (error != CL_SUCCESS)
		return error;

	return CL_SUCCESS;
}

static cl_int _pl_enqueue_projection_kernel(PLContext *pl_ctx, cl_kernel kernel, PLProjection proj, PLProjectionParams *params,
        cl_mem xy_in, cl_mem xy_out, size_t count) {
    cl_int error = CL_SUCCESS;
    int i;
    for (i=0; i<sizeof(_pl_projection_info)/sizeof(_pl_projection_info[0]); i++) {
        if (_pl_projection_info[i].proj == proj) {
            error = _pl_projection_info[i].func(pl_ctx, kernel, params, xy_in, xy_out, count);
            break;
        }
    }
    return error;
}

cl_int pl_enqueue_projection_kernel_points(PLContext *pl_ctx, cl_kernel kernel, PLProjection proj, PLProjectionParams *params,
        PLProjectionBuffer *pl_buf) {
    return _pl_enqueue_projection_kernel(pl_ctx, kernel, proj, params, pl_buf->xy_in, pl_buf->xy_out, pl_buf->count);
}

cl_int pl_enqueue_projection_kernel_grid(PLContext *pl_ctx, cl_kernel kernel, PLProjection proj, PLProjectionParams *params,
        PLPointGridBuffer *src, PLPointGridBuffer *dst) {
    return _pl_enqueue_projection_kernel(pl_ctx, kernel, proj, params, src->grid, dst->grid, src->width * src->height);
}

cl_int pl_enqueue_kernel_albers_equal_area(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
	cl_int error = CL_SUCCESS;
	cl_int argc = 0;
	size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
	PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
	error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &argc);

	double phi1 = params->rlat1 * DEG_TO_RAD;
	double phi2 = params->rlat2 * DEG_TO_RAD;
	double phi0 = params->lat0 * DEG_TO_RAD;

	double rho0;
	double c, n;
    double sinphi, cosphi;

	n = sinphi = sin(phi1);
	cosphi = cos(phi1);

	if (_pl_spheroid_is_spherical(params->spheroid)) {
        n = .5 * (sinphi + sin(phi2));
        c = 1.0 + sin(phi2) * sinphi;
        rho0 = sqrt(c - 2.f * n * sin(phi0));
	} else {
		double ml1, m1;

		m1 = _pl_msfn(sinphi, cosphi, info.ecc2);
		ml1 = _pl_qsfn(sinphi, info.ecc, info.one_ecc2);
		if (fabs(phi1 - phi2) >= EPS7) {
			double ml2, m2;

			sinphi = sin(phi2);
			cosphi = cos(phi2);
			m2 = _pl_msfn(sinphi, cosphi, info.ecc2);
			ml2 = _pl_qsfn(sinphi, info.ecc, info.one_ecc2);
			n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
		}

		c = m1 * m1 + ml1 * n;
		rho0 = sqrt(c - n * _pl_qsfn(sin(phi0), info.ecc, info.one_ecc2));
	}

	if (!_pl_spheroid_is_spherical(params->spheroid)) {
		error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, info.ec);
	}
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->scale * info.major_axis / n);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->y0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->lon0 * DEG_TO_RAD);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, rho0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, c);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, n);

	if (error != CL_SUCCESS)
		return error;

	return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_american_polyconic(PLContext *pl_ctx,cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
	cl_int error;
	cl_int offset = 0;
	size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
	PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
	error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &offset);

	double phi0 = params->lat0 * DEG_TO_RAD;

	double ml0 = _pl_mlfn(phi0, sin(phi0), cos(phi0), info.en);

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->scale * info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->y0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, phi0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->lon0 * DEG_TO_RAD);
    if (!_pl_spheroid_is_spherical(params->spheroid)) {
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, ml0);
        error |= pl_set_kernel_arg_float8(pl_ctx, kernel, offset++, info.en);
    }
	if (error != CL_SUCCESS) {
		return error;
	}

	return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_lambert_azimuthal_equal_area(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
    cl_int error;
    cl_int offset = 0;
    size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
	PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
	error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &offset);

	double phi0 = params->lat0 * DEG_TO_RAD;

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->scale * info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->y0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, phi0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->lon0 * DEG_TO_RAD);
    if (_pl_spheroid_is_spherical(params->spheroid)) {
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, sin(phi0));
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, cos(phi0));
    } else {
        double qp = _pl_qsfn(1.0, info.ecc, info.one_ecc2);
        double sinPhi = sin(phi0);
        double sinB1 = _pl_qsfn(sinPhi, info.ecc, info.one_ecc2) / qp;
        double cosB1 = sqrt(1.0 - sinB1 * sinB1);

        double rq = sqrt(0.5 * qp);
        double dd = cos(phi0) / (sqrt(1.0 - info.ecc2 * sinPhi * sinPhi) * rq * cosB1);
        double ymf = rq / dd;
        double xmf = rq * dd;

        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, qp);
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, sinB1);
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, cosB1);

        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, rq);
        error |= pl_set_kernel_arg_float4(pl_ctx, kernel, offset++, info.apa);

        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, dd);
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, xmf);
        error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, ymf);
    }
    if (error != CL_SUCCESS) {
		return error;
	}

    return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_lambert_conformal_conic(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
    cl_int error;
    cl_int offset = 0;

    size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
    PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
    error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &offset);

    double phi0 = params->lat0 * DEG_TO_RAD;
    double phi1 = params->rlat1 * DEG_TO_RAD;
    double phi2 = params->rlat2 * DEG_TO_RAD;

    double rho0, c, n, sinphi1, cosphi1, sinphi2;
    int secant = 0;

    sinphi1 = sin(phi1);
    cosphi1 = cos(phi1);
    if (fabs(phi1 - phi2) < 1.e-7) {
        n = sinphi1;
    } else {
        secant = 1;
    }

    if (_pl_spheroid_is_spherical(params->spheroid)) {
        if (secant)
            n = log(cosphi1 / cos(phi2)) / (asinh(tan(phi2)) - asinh(tan(phi1)));
        c = cosphi1 * pow(tan(M_PI_4 + .5 * phi1), n) / n;
        rho0 = c * pow(tan(M_PI_4 + .5 * phi0), -n);
    } else {
        double m1, ml1;

        m1 = _pl_msfn(sinphi1, cosphi1, info.ecc2);
        ml1 = _pl_tsfn(phi1, sinphi1, info.ecc);
        if (secant) {
            sinphi2 = sin(phi2);
            n = log(m1 / _pl_msfn(sinphi2, cos(phi2), info.ecc2));
            n /= log(ml1 / _pl_tsfn(phi2, sinphi2, info.ecc));
        }
        c = m1 * pow(ml1, -n) / n;
        rho0 = c * pow(_pl_tsfn(phi0, sin(phi0), info.ecc), n);
    }

    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->scale * info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->y0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->lon0 * DEG_TO_RAD);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, rho0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, c);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, n);

    if (error != CL_SUCCESS) {
		return error;
	}

    return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_mercator(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
	cl_int error;
	cl_int offset = 0;
	size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

	PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
	error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &offset);

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->scale * info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, offset++, params->y0);
	if (error != CL_SUCCESS) {
		return error;
	}
	return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_oblique_stereographic(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
	cl_int error = CL_SUCCESS;
	cl_int argc = 0;
	size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;
	PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);

	error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &argc);

    double phi0 = params->lat0 * DEG_TO_RAD;

    double sinPhi0, cosPhi0;
    sinPhi0 = sin(phi0);
    cosPhi0 = cos(phi0);

    double sinPhiC0, cosPhiC0;
    double scale_r2 = 2.0 * params->scale * info.major_axis * sqrt(info.one_ecc2) / (1.0 - info.ecc2 * sinPhi0 * sinPhi0);

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, scale_r2);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->x0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->y0);

	if (!_pl_spheroid_is_spherical(params->spheroid)) {
        double c0 = sqrt(1.0 + info.ecc2 * cosPhi0 * cosPhi0 * cosPhi0 * cosPhi0 / info.one_ecc2);
        double phiC0 = asin(sinPhi0 / c0);
        sinPhiC0 = sin(phiC0);
        cosPhiC0 = cos(phiC0);

        double k0 = tan(0.5 * phiC0 + M_PI_4) / (
                pow(tan(0.5 * phi0 + M_PI_4), c0) *
                pow((1.-info.ecc * sinPhi0)/(1.+info.ecc*sinPhi0), 0.5 * c0 * info.ecc) );

		error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, c0);
		error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, log(k0));
	} else {
        sinPhiC0 = sinPhi0;
        cosPhiC0 = cosPhi0;
    }

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->lon0 * DEG_TO_RAD);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, sinPhiC0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, cosPhiC0);
	if (error != CL_SUCCESS)
		return error;

	return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_robinson(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
	cl_int error = CL_SUCCESS;
	cl_int argc = 0;

	PLSpheroidInfo info = _pl_get_spheroid_info(PL_SPHEROID_SPHERE);

	error |= pl_set_kernel_arg_mem(pl_ctx, kernel, argc++, xy_in);
	error |= pl_set_kernel_arg_mem(pl_ctx, kernel, argc++, xy_out);
	error |= pl_set_kernel_arg_uint(pl_ctx, kernel, argc++, count);

	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->scale * info.major_axis);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->x0);
	error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->y0);
	if (error != CL_SUCCESS)
		return error;

	return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, count);
}

cl_int pl_enqueue_kernel_transverse_mercator(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
    cl_int error = CL_SUCCESS;
    cl_int argc = 0;
    size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    PLSpheroidInfo info = _pl_get_spheroid_info(params->spheroid);
    error = _pl_set_projection_kernel_args(pl_ctx, kernel, xy_in, xy_out, count, &info, &argc);

    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->scale * info.major_axis * info.krueger_A);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->y0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->lon0 * DEG_TO_RAD);
    if (!_pl_spheroid_is_spherical(params->spheroid)) {
        error |= pl_set_kernel_arg_float8(pl_ctx, kernel, argc++, info.krueger_alpha);
        error |= pl_set_kernel_arg_float8(pl_ctx, kernel, argc++, info.krueger_beta);
    }
    if (error != CL_SUCCESS)
        return error;

    return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_enqueue_kernel_winkel_tripel(PLContext *pl_ctx, cl_kernel kernel,
        PLProjectionParams *params, cl_mem xy_in, cl_mem xy_out, size_t count) {
    cl_int error = CL_SUCCESS;
    int argc = 0;

    size_t vec_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    PLSpheroidInfo info = _pl_get_spheroid_info(PL_SPHEROID_SPHERE);

    double cosphi1 = isnan(params->rlat1) ? M_2_PI : cos(params->rlat1 * DEG_TO_RAD);

    error |= pl_set_kernel_arg_mem(pl_ctx, kernel, argc++, xy_in);
	error |= pl_set_kernel_arg_mem(pl_ctx, kernel, argc++, xy_out);
	error |= pl_set_kernel_arg_uint(pl_ctx, kernel, argc++, vec_count);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->scale * info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->x0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->y0);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, params->lon0 * DEG_TO_RAD);
    error |= pl_set_kernel_arg_float(pl_ctx, kernel, argc++, cosphi1);

    if (error != CL_SUCCESS)
        return error;

    return _pl_enqueue_kernel_1d(pl_ctx->queue, kernel, vec_count);
}

cl_int pl_run_kernel_inverse_geodesic( PLContext *pl_ctx, cl_kernel inv_kernel,
        PLInverseGeodesicBuffer *pl_buf, float *dist_out, PLSpheroid pl_ell, double scale) {
	int argc = 0;
	cl_int error = CL_SUCCESS;
	PLSpheroidInfo info = _pl_get_spheroid_info(pl_ell);

	size_t xy2VecCount = ck_padding(pl_buf->xy2_count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

	error |= pl_set_kernel_arg_mem(pl_ctx, inv_kernel, argc++, pl_buf->xy1_in);
	error |= pl_set_kernel_arg_mem(pl_ctx, inv_kernel, argc++, pl_buf->xy2_in);
	error |= pl_set_kernel_arg_mem(pl_ctx, inv_kernel, argc++, pl_buf->dist_out);
	error |= pl_set_kernel_arg_float(pl_ctx, inv_kernel, argc++, info.major_axis);

	if (error != CL_SUCCESS) {
		return error;
	}

	const size_t dim[2] = { pl_buf->xy1_count, xy2VecCount };

	error = clEnqueueNDRangeKernel(pl_ctx->queue, inv_kernel, 2, NULL, dim, NULL, 0, NULL, NULL);
	if (error != CL_SUCCESS) {
		return error;
	}

	float *dist_pad_out;
	size_t pad_count = 2 * pl_buf->xy1_count * ck_padding(pl_buf->xy2_count, PL_FLOAT_VECTOR_SIZE);
	if ((dist_pad_out = malloc(sizeof(float) * pad_count)) == NULL) {
		return CL_OUT_OF_HOST_MEMORY;
	}

	error = clEnqueueReadBuffer(pl_ctx->queue, pl_buf->dist_out, CL_TRUE, 0, pad_count * sizeof(cl_float),
								dist_pad_out, 0, NULL, NULL);
	if (error != CL_SUCCESS) {
		free(dist_pad_out);
		return error;
	}

	error = clFinish(pl_ctx->queue);
	if (error != CL_SUCCESS) {
		free(dist_pad_out);
		return error;
	}

	int j_size = ck_padding(pl_buf->xy2_count, PL_FLOAT_VECTOR_SIZE);
    int i, j;

	for (i=0; i<pl_buf->xy1_count; i++) {
		for (j=0; j<pl_buf->xy2_count; j++) {
			dist_out[i*pl_buf->xy2_count+j] = dist_pad_out[i*j_size+j] * scale;
		}
	}

	free(dist_pad_out);

	return CL_SUCCESS;
}

cl_int pl_run_kernel_forward_geodesic_fixed_distance(PLContext *pl_ctx, cl_kernel fwd_kernel,
    PLForwardGeodesicFixedDistanceBuffer *pl_buf, float *xy_out, PLSpheroid pl_ell, double distance) {
	PLSpheroidInfo info = _pl_get_spheroid_info(pl_ell);

	int argc = 0;

	cl_int error = CL_SUCCESS;

	double D, sinD, cosD;

	D = distance / info.major_axis;
	sinD = sin(D);
	cosD = cos(D);

	error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->xy_in);
	error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->phi_sincos);
	error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->az_sincos);
	error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->xy_out);
	error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, D);
	error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, sinD);
	error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, cosD);
	if (!_pl_spheroid_is_spherical(pl_ell)) {
		double flattening = 1.f/info.inverse_flattening;
		error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, flattening);
	}
	if (error != CL_SUCCESS) {
		return error;
	}

	const size_t dim[2] = { pl_buf->xy_count,
		ck_padding(pl_buf->az_count, PL_FLOAT_VECTOR_SIZE)/PL_FLOAT_VECTOR_SIZE };

	error = clEnqueueNDRangeKernel(pl_ctx->queue, fwd_kernel, 2, NULL,
								   dim, NULL, 0, NULL, NULL);
	if (error != CL_SUCCESS) {
		return error;
	}

	error = clEnqueueReadBuffer(pl_ctx->queue, pl_buf->xy_out, CL_TRUE, 0,
								2 * pl_buf->xy_count * pl_buf->az_count * sizeof(cl_float),
								xy_out, 0, NULL, NULL);
	if (error != CL_SUCCESS) {
		return error;
	}

	error = clFinish(pl_ctx->queue);
	if (error != CL_SUCCESS) {
		return error;
	}

	return CL_SUCCESS;
}

cl_int pl_run_kernel_forward_geodesic_fixed_angle(PLContext *pl_ctx, cl_kernel fwd_kernel,
    PLForwardGeodesicFixedAngleBuffer *pl_buf, double xy_in[2], float *xy_out, PLSpheroid pl_ell, double angle) {
    int argc = 0;

    size_t distVecCount = ck_padding(pl_buf->dist_count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    cl_int error = CL_SUCCESS;

    double sinAzimuth = sin(angle);
    double cosAzimuth = cos(angle);

    error |= pl_set_kernel_arg_float2(pl_ctx, fwd_kernel, argc++, xy_in);
    error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->dist_in);
    error |= pl_set_kernel_arg_mem(pl_ctx, fwd_kernel, argc++, pl_buf->xy_out);
    error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, angle);
    error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, sinAzimuth);
    error |= pl_set_kernel_arg_float(pl_ctx, fwd_kernel, argc++, cosAzimuth);
    if (error != CL_SUCCESS) {
        return error;
    }

    error = clEnqueueNDRangeKernel(pl_ctx->queue, fwd_kernel, 1, NULL, &distVecCount,
                                   NULL, 0, NULL, NULL);
    if (error != CL_SUCCESS) {
        return error;
    }

    error = clEnqueueReadBuffer(pl_ctx->queue, pl_buf->xy_out, CL_TRUE, 0,
                                2 * pl_buf->dist_count * sizeof(cl_float),
                                xy_out, 0, NULL, NULL);
    if (error != CL_SUCCESS) {
        return error;
    }

    error = clFinish(pl_ctx->queue);
    if (error != CL_SUCCESS) {
        return error;
    }

    return CL_SUCCESS;
}

cl_int pl_run_kernel_geodesic_to_cartesian(PLContext *pl_ctx, cl_kernel g2c_kernel,
        PLDatumShiftBuffer *pl_buf, PLSpheroid pl_ell) {
    PLSpheroidInfo info = _pl_get_spheroid_info(pl_ell);

    int argc = 0;
    cl_int error = CL_SUCCESS;
    size_t vec_count = ck_padding(pl_buf->count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    error |= pl_set_kernel_arg_mem(pl_ctx, g2c_kernel, argc++, pl_buf->xy_in);
    error |= pl_set_kernel_arg_mem(pl_ctx, g2c_kernel, argc++, pl_buf->x_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, g2c_kernel, argc++, pl_buf->y_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, g2c_kernel, argc++, pl_buf->z_rw);

    error |= pl_set_kernel_arg_float(pl_ctx, g2c_kernel, argc++, info.ecc);
    error |= pl_set_kernel_arg_float(pl_ctx, g2c_kernel, argc++, info.ecc2);
    error |= pl_set_kernel_arg_float(pl_ctx, g2c_kernel, argc++, info.one_ecc2);

    error |= pl_set_kernel_arg_float(pl_ctx, g2c_kernel, argc++, info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, g2c_kernel, argc++, info.minor_axis);

    if (error != CL_SUCCESS)
        return error;

    error = clEnqueueNDRangeKernel(pl_ctx->queue, g2c_kernel, 1, NULL, &vec_count, NULL, 0, NULL, NULL);
    if (error != CL_SUCCESS)
        return error;

    error = clFinish(pl_ctx->queue);
    if (error != CL_SUCCESS)
        return error;

    return CL_SUCCESS;
}

/* We cut the computations in half by doing some matrix algebra beforehand.
 * Each datum shift is essentially an affine transformation in 3D cartesian
 * coordinates. So to get from on datum to another, instead of transforming
 * to/from WGS 84, we concatenate the transformation matrix of the source
 * datum to the inverse transformation matrix of the destination matrix, and
 * then just do a single matrix multiplication on each point instead of two
 * in order to transform it.
 */
cl_int pl_run_kernel_transform_cartesian(PLContext *pl_ctx, cl_kernel transform_kernel,
        PLDatumShiftBuffer *pl_buf, PLDatum src_datum, PLDatum dst_datum) {
    double Rx1, Ry1, Rz1, M1, Dx1, Dy1, Dz1;
    double Rx2, Ry2, Rz2, M2, Dx2, Dy2, Dz2;

    cl_int error = CL_SUCCESS;
    int argc = 0;
    size_t vec_count = ck_padding(pl_buf->count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    Dx1 = pl_datum_params[src_datum].dx;
    Dy1 = pl_datum_params[src_datum].dy;
    Dz1 = pl_datum_params[src_datum].dz;
    M1  = pl_datum_params[src_datum].ppm*1.e-6 + 1;
    Rx1 = pl_datum_params[src_datum].ex * SEC_TO_RAD;
    Ry1 = pl_datum_params[src_datum].ey * SEC_TO_RAD;
    Rz1 = pl_datum_params[src_datum].ez * SEC_TO_RAD;

    Dx2 = pl_datum_params[dst_datum].dx;
    Dy2 = pl_datum_params[dst_datum].dy;
    Dz2 = pl_datum_params[dst_datum].dz;
    M2  = pl_datum_params[dst_datum].ppm*1.e-6 + 1;
    Rx2 = pl_datum_params[dst_datum].ex * SEC_TO_RAD;
    Ry2 = pl_datum_params[dst_datum].ey * SEC_TO_RAD;
    Rz2 = pl_datum_params[dst_datum].ez * SEC_TO_RAD;

    struct pl_matrix matrix1 = pl_affine_transform_make(Rx1, Ry1, Rz1, M1, Dx1, Dy1, Dz1);
    struct pl_matrix matrix2 = pl_affine_transform_make(Rx2, Ry2, Rz2, M2, Dx2, Dy2, Dz2);

    /* Invert the destination matrix */
    __CLPK_integer n = 4;
    __CLPK_integer info;

    __CLPK_integer ipiv[4];
    __CLPK_doublereal work[4];
    __CLPK_integer work_n = 4;

    dgetrf_(&n, &n, &matrix2.elements[0][0], &n, ipiv, &info);
    if (info != 0) {
        return info;
    }

    dgetri_(&n, &matrix2.elements[0][0], &n, ipiv, work, &work_n, &info);
    if (info != 0) {
        return info;
    }

    __CLPK_doublereal alpha = 1.0;
    __CLPK_doublereal beta = 0.0;
    __CLPK_doublereal result_matrix[4][4];

    /* Multiply the source matrix with the inverse destination matrix */
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha,
          &matrix2.elements[0][0], n,
          &matrix1.elements[0][0], n,
          beta, &result_matrix[0][0], n);

    /* transpose the result and store single-precision */
    double tmatrix[4][4];
    int i, j;
    for (i=0; i<4; i++) {
        for (j=0; j<4; j++) {
            tmatrix[i][j] = result_matrix[j][i];
        }
    }

    error |= pl_set_kernel_arg_mem(pl_ctx, transform_kernel, argc++, pl_buf->x_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, transform_kernel, argc++, pl_buf->y_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, transform_kernel, argc++, pl_buf->z_rw);
    error |= pl_set_kernel_arg_float16(pl_ctx, transform_kernel, argc++, &tmatrix[0][0]);
    if (error != CL_SUCCESS)
        return error;

    error = clEnqueueNDRangeKernel(pl_ctx->queue, transform_kernel, 1, NULL,
                                   &vec_count, NULL, 0, NULL, NULL);
    if (error != CL_SUCCESS)
        return error;

    return error;
}

cl_int pl_run_kernel_cartesian_to_geodesic(PLContext *pl_ctx, cl_kernel c2g_kernel,
        PLDatumShiftBuffer *pl_buf, float *xy_out, PLSpheroid pl_ell) {
    PLSpheroidInfo info = _pl_get_spheroid_info(pl_ell);

    int argc = 0;
    cl_int error = CL_SUCCESS;
    size_t vec_count = ck_padding(pl_buf->count, PL_FLOAT_VECTOR_SIZE) / PL_FLOAT_VECTOR_SIZE;

    error |= pl_set_kernel_arg_mem(pl_ctx, c2g_kernel, argc++, pl_buf->x_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, c2g_kernel, argc++, pl_buf->y_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, c2g_kernel, argc++, pl_buf->z_rw);
    error |= pl_set_kernel_arg_mem(pl_ctx, c2g_kernel, argc++, pl_buf->xy_out);

    error |= pl_set_kernel_arg_float(pl_ctx, c2g_kernel, argc++, info.ecc);
    error |= pl_set_kernel_arg_float(pl_ctx, c2g_kernel, argc++, info.ecc2);
    error |= pl_set_kernel_arg_float(pl_ctx, c2g_kernel, argc++, info.one_ecc2);

    error |= pl_set_kernel_arg_float(pl_ctx, c2g_kernel, argc++, info.major_axis);
    error |= pl_set_kernel_arg_float(pl_ctx, c2g_kernel, argc++, info.minor_axis);
    if (error != CL_SUCCESS)
        return error;

    error = clEnqueueNDRangeKernel(pl_ctx->queue, c2g_kernel, 1, NULL,
                                   &vec_count, NULL, 0, NULL, NULL);
    if (error != CL_SUCCESS)
        return error;

    if (xy_out != NULL) {
        error = clEnqueueReadBuffer(pl_ctx->queue, pl_buf->xy_out, CL_TRUE, 0,
                                    2 * pl_buf->count * sizeof(cl_float), xy_out, 0, NULL, NULL);
        if (error != CL_SUCCESS)
            return error;
    }

    error = clFinish(pl_ctx->queue);
    if (error != CL_SUCCESS)
        return error;

    return CL_SUCCESS;
}

ext: c

max_stars_repo_name: evanmiller/ProjCL

max_issues_repo_name: evanmiller/ProjCL

max_forks_repo_name: evanmiller/ProjCL


text:
/* Internal header file for cminpack, by Frederic Devernay. */
#ifndef __CMINPACKP_H__
#define __CMINPACKP_H__

#ifndef __CMINPACK_H__
#error "cminpackP.h in an internal cminpack header, and must be included after all other headers (including cminpack.h)"
#endif

#if (defined (USE_CBLAS) || defined (USE_LAPACK)) && !defined (__cminpack_double__)
#error "cminpack can use cblas and lapack only in double precision mode"
#endif

#ifdef USE_CBLAS
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif
#define __cminpack_enorm__(n,x) cblas_dnrm2(n,x,1)
#else
#define __cminpack_enorm__(n,x) __cminpack_func__(enorm)(n,x)
#endif

#ifdef USE_LAPACK
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
#if defined(__LP64__) /* In LP64 match sizes with the 32 bit ABI */
typedef int 		__CLPK_integer;
typedef int 		__CLPK_logical;
typedef float 		__CLPK_real;
typedef double 		__CLPK_doublereal;
typedef __CLPK_logical 	(*__CLPK_L_fp)();
typedef int 		__CLPK_ftnlen;
#else
typedef long int 	__CLPK_integer;
typedef long int 	__CLPK_logical;
typedef float 		__CLPK_real;
typedef double 		__CLPK_doublereal;
typedef __CLPK_logical 	(*__CLPK_L_fp)();
typedef long int 	__CLPK_ftnlen;
#endif
//extern void dlartg_(double *f, double *g, double *cs, double *sn, double *r__);
int dlartg_(__CLPK_doublereal *f, __CLPK_doublereal *g, __CLPK_doublereal *cs,
            __CLPK_doublereal *sn, __CLPK_doublereal *r__)
//extern void dgeqp3_(int *m, int *n, double *a, int *lda, int *jpvt, double *tau, double *work, int *lwork, int *info);
int dgeqp3_(__CLPK_integer *m, __CLPK_integer *n, __CLPK_doublereal *a, __CLPK_integer *
            lda, __CLPK_integer *jpvt, __CLPK_doublereal *tau, __CLPK_doublereal *work, __CLPK_integer *lwork,
            __CLPK_integer *info)
//extern void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, double *work, int *lwork, int *info);
int dgeqrf_(__CLPK_integer *m, __CLPK_integer *n, __CLPK_doublereal *a, __CLPK_integer *
            lda, __CLPK_doublereal *tau, __CLPK_doublereal *work, __CLPK_integer *lwork, __CLPK_integer *info)
#endif
#endif

#define real_mp __cminpack_real__
#define minpack_min(a,b) ((a) <= (b) ? (a) : (b))
#define minpack_max(a,b) ((a) >= (b) ? (a) : (b))
#define TRUE_ (1)
#define FALSE_ (0)

#endif /* !__CMINPACKP_H__ */

ext: h

max_stars_repo_name: ptrbortolotti/WISDEM

max_issues_repo_name: ptrbortolotti/WISDEM

max_forks_repo_name: ptrbortolotti/WISDEM


text:
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_chebyshev.h>

double
f (double x, void *p)
{
  p = 0;
  if (x < 0.5)
    return 0.25;
  else
    return 0.75;
}

int
main (void)
{
  double x;

  gsl_cheb_series *cs = gsl_cheb_alloc (40);

  gsl_function F;

  F.function = f;
  F.params = 0;

  gsl_cheb_init (cs, &F, 0.0, 1.0);

  for (x = -1; x < 2; x += 0.001)
    {
      double r10 = gsl_cheb_eval_n (cs, 10, x);
      double r40 = gsl_cheb_eval (cs, x);
      printf ("%g %g %g %g\n", x, GSL_FN_EVAL (&F, x), r5, r40);
    }

  gsl_cheb_free (cs);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <petsc.h>
#include <petsc/private/fortranimpl.h>
#include <omp.h>
//#include <papi.h>
/*
  Global Variables
*/
double  gradTime; // Gradient calculation time
double  fluxTime; // Flux evaluation timing
int     nThreads; // Number of threads at each MPI process
int     rank;     // MPI process ID
int     size;     // MPI number of running processes

#ifdef HW_COUNTER
int       eventset; // PAPI counter context
long long handler;  // PAPI counter return value
#endif

/*
  FORTRAN routine definitions
*/
#if defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
#define f77name(ucase,lcase,lcbar) lcbar
#elif defined(PETSC_HAVE_FORTRAN_CAPS)
#define f77name(ucase,lcase,lcbar) ucase
#else
#define f77name(ucase,lcase,lcbar) lcase
#endif
/*
  Node coordinates in a 3D mesh
*/
struct COORD
{
  double *x;  // X component
  double *y;  // Y component
  double *z;  // Z component
};
/*
  Normal to faces and length
*/
struct NORMAL
{
  double  *len;         // Length
  struct  COORD coord;  // X, Y, and Z components
};
/*
  Weights for calculating the gradients using least square
*/
struct WEIGHT
{
  double *r11;  // w11
  double *r12;  // r12r11
  double *r13;  // r13r11
  double *r22;  // w22
  double *r23;  // r23r22
  double *r33;  // w33
  double *r44;  // rmess
};
/*
  MeTiS partitioning using for OpenMP
*/
struct THREAD
{
  int     nedge;          // Number edge for all threads with replication
  int     *partitions;    // MeTiS nodes partitioning
  int     *nedgeLoc;      // Number of edges local for each thread
  int     *eptr;          // Edge pointers
  struct  NORMAL  normal; // Normal to faces and length
};
/*
  A boundary's elements
*/
struct BOUNDELM
{
  int     nnode;        // Number of nodes
  int     nfacet;       // Number of facets
  int     *f2nt;        // Facet to tetrahedron indices
  int     *inode;       // nodes pointers
  struct  COORD normal; // Normal to faces
};
/*
  Boundaries
*/
struct BOUND
{
  int     nsbound;            // Number of solid boundaries
  int     nvbound;            // Number of viscous boundaries
  int     nfbound;            // Number of free boundaries
  struct  BOUNDELM solid;     // Solid boundary elements
  struct  BOUNDELM viscs;     // Viscous boundary elements
  struct  BOUNDELM ffild;     // Free boundary elements
};
/*
  Grid information
*/
typedef struct grid_t
{
  int     nnodes;           // Number of mesh nodes
  int     nedge;            // Number of mesh edges
  int     *eptr;            // Edge pointers array
  struct  NORMAL  normal;   // Normal to faces and length
  struct  COORD   coord;    // Node Coordinates
  struct  THREAD  thread;   // MeTiS subdomain partitioning elements
  struct  BOUND   bound;    // Boundaries elements
  struct  WEIGHT  weight;   // Least square weights
  double  *area;            // Area of control volumes
  double  *cdt;             // Pseudo-time step for each cell
} GRID;
/*
  PETSc data structures
*/
typedef struct gridPETSc_t
{
  Vec   qnode;        // Global distributed solution vector
  Vec   qnodeLoc;     // Local sequential solution vector
  Vec   res;          // Residual vector
  Vec   grad;         // Gradient vector
  Mat   A;            // Left hand side matrix
} GRIDPETSc;
/*
  Time step context
*/
typedef struct tstepctx_t
{
  /*
    A vector state at which to evaluate residual (x) at every
    time step (hold the old x vector before proceeding to the
    next SNES time step)
  */
  Vec     qnode;
  /*
    A vector to put residual (function value) (f) used at every
    time step calculation
  */
  Vec     res;
  /*
    Courant–Friedrichs–Lewy stability restriction condition
  */
  double  CFL;
  double  CFLInit;
  double  CFLMax;
  /*
    2-norm of function at current iterate
  */
  double  fnorm;
  double  fnormInit;
  double  fnormRatio;
  /*
    A flag to indicate that the formFunction is called to compute the
    pseudo time step
  */
  int     isPseudoTimeStep;
  /*
    Time step counting variable: time step loop iterator (i)
  */
  int     iTimeStep;
  /*
    Maximum time steps for the nonlinear solve
  */
  int     maxTimeSteps;
  /*
    A flag for local or global time step
  */
  int     isLocal;
} TstepCtx;
/*
  Encapsulate the grid data structure information and time step
  data structure information
*/
typedef struct appctx_t
{
  GRID        *grid;        // Pointer to grid info
  GRIDPETSc   *gridPETSc;   // PETSc data structure context
  TstepCtx    *tsCtx;       // Pointer to time stepping context
  PetscBool   isPreLoading; // Preloading flag (single time step)
} AppCtx;


text:
#ifndef KGS_GSLSVD_H
#define KGS_GSLSVD_H

#include "math/SVD.h"
#include <gsl/gsl_matrix.h>

class SVDGSL: public SVD {
 public:
  SVDGSL(gsl_matrix* M): SVD(M){}

 protected:

  void UpdateFromMatrix() override;
};

#endif //KGS_GSLSVD_H


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <pthread.h>

#include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

#include "LNAbundance.h"

static int  verbose  = FALSE;

static int  nLines   = 9;

static char *usage[] = {"LNAbundance - \n",
                        "Required parameters:\n",
			"   -a    abundance data\n",
			"   -in   filename      parameter file  \n",
                        "Optional:\n",
			"   -b    integer       length of sample file to ignore\n",
			"   -s    integer       sampling frequency\n",
			"   -seed long          seed random number generator\n",
			"   -v                  verbose\n"};

void updateExpectations(double* adExpect, int nMax, double dMDash, double dV, int nS, t_Data *ptData)
{
  int i = 0;

  for(i = 1; i <= nMax; i++){
    int nA = i;
    double dLog = 0.0, dP = 0.0;

    if(nA < MAX_QUAD){
      dLog = logLikelihoodQuad(nA, dMDash, dV);
    }
    else{
      dLog = logLikelihoodRampal(nA, dMDash, dV);
    }

    dP = exp(dLog);

    adExpect[i - 1]+= dP*nS;
  }
}

int main(int argc, char* argv[]){
  t_Params    tParams;
  t_LNParams* atLNParams;
  int         i = 0, nSamples = 0, nMax = 0;
  t_Data      tData;
  double*     adExpect = NULL;

  gsl_set_error_handler_off();

  getCommandLineParams(&tParams, argc, argv);

  /*allocate memory for samples*/
  atLNParams = (t_LNParams *) malloc(MAX_SAMPLES*sizeof(t_LNParams));
  if(!atLNParams)
    goto memoryError;

  /*read in Monte-Carlo samples*/
  readSamples(&tParams, atLNParams, &nSamples);

  readAbundanceData(tParams.szAbundFile, &tData);

  nMax = tData.aanAbund[tData.nNA - 1][0];
  nMax = floor(pow(2.0,ceil(log((double) nMax)/log(2.0)) + 2.0) + 1.0e-7);

  adExpect = (double *) malloc(sizeof(double)*nMax);

  for(i = 0; i < nMax; i++){
    adExpect[i] = 0.0;
  }

  for(i = 0; i < nSamples; i++){
    updateExpectations(adExpect, nMax,
		       atLNParams[i].dMDash,
		       atLNParams[i].dV,
		       atLNParams[i].nS,
		       &tData);
  }

  for(i = 1; i <= nMax; i++){
    printf("%d %f\n",i,adExpect[i - 1]/((double) nSamples));
  }

  free(adExpect);
  exit(EXIT_SUCCESS);

 memoryError:
  fprintf(stderr, "Failed to allocate memory in main aborting ...\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int intCompare(const void *pvA, const void *pvB)
{
  int* pnA = (int *) pvA, *pnB = (int *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

int doubleCompare(const void *pvA, const void *pvB)
{
  double* pnA = (double *) pvA, *pnB = (double *) pvB;

  if(*pnA < *pnB)
    return +1;
  else if(*pnA == *pnB){
    return 0;
  }
  else{
    return -1;
  }
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

 /*get abundance file name*/
  ptParams->szAbundFile  = extractParameter(argc,argv, ABUND_FILE,ALWAYS);
  if(ptParams->szAbundFile == NULL)
    goto error;

  /*get long seed*/
  szTemp  = extractParameter(argc,argv,SEED,OPTION);
  if(szTemp != NULL){
    ptParams->lSeed = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->lSeed = 0;
  }

  /*get burn*/
  szTemp  = extractParameter(argc, argv, BURN, OPTION);
  if(szTemp != NULL){
    ptParams->nBurn = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nBurn = DEF_BURN;
  }

 /*get long seed*/
  szTemp  = extractParameter(argc, argv, SAMPLE, OPTION);
  if(szTemp != NULL){
    ptParams->nSample = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nSample = DEF_SAMPLE;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readSamples(t_Params *ptParams, t_LNParams *atLNParams, int *pnSamples)
{
  int  nSamples = 0;
  char *szInputFile = ptParams->szInputFile;
  char szLine[MAX_LINE_LENGTH];
  FILE *ifp = NULL;

  ifp = fopen(szInputFile, "r");

  if(ifp){
    while(fgets(szLine, MAX_LINE_LENGTH, ifp)){
      char *szTok = NULL, *szBrk = NULL, *pcError = NULL;
      int  nTime  = 0;

      /*remove trailing new line*/
      szBrk = strpbrk(szLine, "\n"); (*szBrk) = '\0';

      szTok = strtok(szLine, DELIM);

      nTime = strtol(szTok, &pcError, 10);
      if(*pcError != '\0') goto fileFormatError;

      if(nTime > ptParams->nBurn && nTime % ptParams->nSample == 0){

	szTok = strtok(NULL,DELIM);
	atLNParams[nSamples].dMDash = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atLNParams[nSamples].dV = strtod(szTok, &pcError);
	if(*pcError != '\0') goto fileFormatError;

	szTok = strtok(NULL,DELIM);
	atLNParams[nSamples].nS = strtol(szTok, &pcError, 10);
	if(*pcError != '\0') goto fileFormatError;

	nSamples++;
      }
    }

    fclose(ifp);
  }
  else{
    fprintf(stderr, "Failed to open file %s aborting\n", szInputFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  (*pnSamples) = nSamples;
  return;

 fileFormatError:
  fprintf(stderr, "Incorrectly formatted input file aborting\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

double f1(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
  double dTemp = (x - dMDash);
  double dExp  = x*((double) n) - exp(x) - 0.5*((dTemp*dTemp)/dV);
  double dRet  = exp(dExp);

  return dRet;
}

double logLikelihoodRampal(int n, double dMDash, double dV)
{
  double dN = (double) n;
  double dLogLik = 0.0, dTemp = gsl_pow_int(log(dN) - dMDash,2), dTemp3 = gsl_pow_int(log(dN) - dMDash,3);

  dLogLik = -0.5*log(2.0*M_PI*dV) - log(dN) - (dTemp/(2.0*dV));

  dLogLik += log(1.0 + 1.0/(2.0*dN*dV)*(dTemp/dV + log(dN) - dMDash - 1.0)
		 + 1.0/(6.0*dN*dN*dV*dV*dV)*(3.0*dV*dV - (3.0*dV - 2.0*dV*dV)*(dMDash - log(dN))
		 - 3.0*dV*dTemp + dTemp3));

  return dLogLik;
}

double logLikelihoodQuad(int n, double dMDash, double dV)
{
  gsl_integration_workspace *ptGSLWS =
    gsl_integration_workspace_alloc(1000);
  double dLogFac1   = 0.0, dLogFacN  = 0.0;
  double dResult = 0.0, dError = 0.0, dPrecision = 0.0;
  gsl_function tGSLF;
  t_LNParams tLNParams;
  double dEst = dMDash + ((double)n)*dV, dA = 0.0, dB = 0.0;

  tLNParams.n = n; tLNParams.dMDash = dMDash; tLNParams.dV = dV;

  tGSLF.function = &f1;
  tGSLF.params   = (void *) &tLNParams;

  dLogFac1 = log(2.0*M_PI*dV);

  if(n < 50){
    dLogFacN = gsl_sf_fact(n);
    dLogFacN = log(dLogFacN);
  }
  else{
    dLogFacN = gsl_sf_lngamma(((double) n) + 1.0);
  }

  if(dEst > dV){
    double dMax = 0.0;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    solveF(0.0, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);

    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }
  else{
    double dMax = 0.0;
    double dLower = dEst - dV;
    double dUpper = (((double) n) + (dMDash/dV) - 1.0)/(1.0 + 1.0/dV);
    double dVar   = 0.0;

    solveF(dLower, dUpper, derivExponent, (void *) &tLNParams, 1.0e-5, &dMax);

    dVar = sqrt(1.0/((1.0/dV) + exp(dMax)));

    dA = dMax - V_MULT*dVar; dB = dMax + V_MULT*dVar;
  }

  if(n < 10){
    dPrecision = HI_PRECISION;
  }
  else{
    dPrecision = LO_PRECISION;
  }

  gsl_integration_qag(&tGSLF, dA, dB, dPrecision, 0.0, 1000, GSL_INTEG_GAUSS61, ptGSLWS, &dResult, &dError);

  gsl_integration_workspace_free(ptGSLWS);

  return log(dResult) - dLogFacN -0.5*dLogFac1;
}

double derivExponent(double x, void *pvParams)
{
  t_LNParams *ptLNParams = (t_LNParams *) pvParams;
  double dMDash = ptLNParams->dMDash, dV = ptLNParams->dV, n = ptLNParams->n;
  double dTemp = (x - dMDash)/dV, dRet = 0.0;

  dRet = ((double) n) - exp(x) - dTemp;

  return dRet;
}

int solveF(double x_lo, double x_hi, double (*f)(double, void*),
	   void* params, double tol, double *xsolve)
{
  int status, iter = 0, max_iter = 100;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double r = 0;
  gsl_function F;

  F.function = f;
  F.params = params;

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);

  do{
    iter++;
    status = gsl_root_fsolver_iterate (s);
    r = gsl_root_fsolver_root (s);
    x_lo = gsl_root_fsolver_x_lower (s);
    x_hi = gsl_root_fsolver_x_upper (s);

    status = gsl_root_test_interval (x_lo, x_hi, 0, tol);
  }
  while (status == GSL_CONTINUE && iter < max_iter);

  (*xsolve) = gsl_root_fsolver_root (s);
  gsl_root_fsolver_free (s);

  return status;
}

void readAbundanceData(const char *szFile, t_Data *ptData)
{
  int **aanAbund = NULL;
  int  i = 0, nNA = 0, nA = 0, nC = 0;
  int  nL = 0, nJ = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    szTok = strtok(szLine, DELIM2);

    nNA = strtol(szTok,&pcError,10);
    if(*pcError != '\0'){
      goto formatError;
    }

    aanAbund = (int **) malloc(nNA*sizeof(int*));

    for(i = 0; i < nNA; i++){
      aanAbund[i] = (int *) malloc(sizeof(int)*2);

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM2);

      nA = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      szTok = strtok(NULL, DELIM2);

      nC = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      nL += nC;
      nJ += nC*nA;

      aanAbund[i][0]  = nA;
      aanAbund[i][1]  = nC;
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  ptData->nJ          = nJ;
  ptData->nL          = nL;
  ptData->aanAbund    = aanAbund;
  ptData->nNA         = nNA;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_fit.h>

#include "cosmocalc.h"

double transfer_function(double k)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static gsl_spline *cosmocalc_transfer_function_spline = NULL;
  static gsl_interp_accel *cosmocalc_transfer_function_acc = NULL;
  static double c0,c1;

  double transfer_function_table[COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH];
  double k_table[COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH];
  long i;
  double cov00,cov01,cov11,sumsq;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      for(i=0;i<COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH;++i)
        {
          k_table[i] = log(K_MAX/K_MIN)/(COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH-1.0)*((double) i) + log(K_MIN);
          if(cosmoData.useSmoothTransFunc)
	    transfer_function_table[i] = log(transfunct_eh98_smooth(exp(k_table[i])));
	  else if(cosmoData.useNoTransFunc)
	    transfer_function_table[i] = 0.0;
	  else
	    transfer_function_table[i] = log(transfunct_eh98(exp(k_table[i])));
        }

      //init the spline and accelerators
      if(cosmocalc_transfer_function_spline != NULL)
        gsl_spline_free(cosmocalc_transfer_function_spline);
      cosmocalc_transfer_function_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_transfer_function_spline,k_table,transfer_function_table,(size_t) (COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH));
      if(cosmocalc_transfer_function_acc != NULL)
        gsl_interp_accel_reset(cosmocalc_transfer_function_acc);
      else
        cosmocalc_transfer_function_acc = gsl_interp_accel_alloc();

      gsl_fit_linear(k_table+COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH-COSMOCALC_TRANSFER_FUNCTION_FIT_LENGTH,(size_t) 1,
                     transfer_function_table+COSMOCALC_TRANSFER_FUNCTION_TABLE_LENGTH-COSMOCALC_TRANSFER_FUNCTION_FIT_LENGTH,(size_t) 1,
                     (size_t) COSMOCALC_TRANSFER_FUNCTION_FIT_LENGTH,&c0,&c1,&cov00,&cov01,&cov11,&sumsq);
    }

  if(k < K_MIN)
    return 1.0;
  else if(k < K_MAX)
    return exp(gsl_spline_eval(cosmocalc_transfer_function_spline,log(k),cosmocalc_transfer_function_acc));
  else
    return exp(c0+c1*log(k));
}

#define USEPAPER_EH98

double transfunct_eh98(double kin)
{
  //vars
  double k;
  double Tk,Tb,Tc,T0t,T0t11,T0t1bc;
  double omb,om0,omc,h;
  double theta2p7,f,ac,bc,s,q;
  double C,C11,C1bc,st,bb,ab,ksilk,bnode,Gy,y,zeq;
  double a1,a2,b1,b2,keq,Rd,zd,Req,b1d,b2d;

  //get cosmoparms
  omb = cosmoData.OmegaB;
  om0 = cosmoData.OmegaM;
  omc = cosmoData.OmegaM - cosmoData.OmegaB;
  h = cosmoData.h;

  //convert k from hMpc^-1 to Mpc^-1
  k = kin*h;

  //-----------
  //input parms
  //-----------
  theta2p7 = TCMB/2.7;

  //eqn 2
  zeq = 2.50e4*om0*h*h/(theta2p7*theta2p7*theta2p7*theta2p7);

  //eqn 3
  keq = 7.46e-2*om0*h*h/(theta2p7*theta2p7); //Mpc^-{1} (NOT h/Mpc)

  //eqn 4
  b1d = 0.313*pow(om0*h*h,-0.419)*(1.0 + 0.607*pow(om0*h*h,0.674));
  b2d = 0.238*pow(om0*h*h,0.223);
  zd = 1291.0*pow(om0*h*h,0.251)/(1.0 + 0.659*pow(om0*h*h,0.828))
    *(1.0 + b1d*pow(omb*h*h,b2d));

  //eqn 5
#ifdef USEPAPER_EH98
  Rd = 31.5*omb*h*h
    /(theta2p7*theta2p7*theta2p7*theta2p7)
    /((zd)/1e3); //matches wayne's paper, but not the code (zd -> 1+zd)
  //fprintf(stderr,"paper eh 98 1\n");
#else
  Rd = 31.5*omb*h*h
    /(theta2p7*theta2p7*theta2p7*theta2p7)
    /((1.0+zd)/1e3);
#endif
Req = 31.5*omb*h*h/(theta2p7*theta2p7*theta2p7*theta2p7)/(zeq/1e3);

  //eqn 6
  s = 2.0/3.0/keq*sqrt(6.0/Req)*
    log((sqrt(1.0 + Rd) + sqrt(Rd + Req))/(1.0 + sqrt(Req)));

  //eqn 7
  ksilk = 1.6*pow(omb*h*h,0.52)*pow(om0*h*h,0.73)*(1.0 + pow(10.4*om0*h*h,-0.95));

  //eqn 10
  q = k/13.41/keq;

  //eqn 11
  a1 = pow(46.9*om0*h*h,0.670)*(1.0 + pow(32.1*om0*h*h,-0.532));
  a2 = pow(12.0*om0*h*h,0.424)*(1.0 + pow(45.0*om0*h*h,-0.582));
  ac = pow(a1,-1.0*omb/om0)*pow(a2,-1.0*(omb/om0)*(omb/om0)*(omb/om0));

  //eqn 12
  b1 = 0.944/(1.0 + pow(458.0*om0*h*h,-0.708));
  b2 = pow(0.395*om0*h*h,-0.0266);
  bc = 1.0/(1.0 + b1*(pow(omc/om0,b2) - 1.0));

  //eqn 15
#ifdef USEPAPER_EH98
  y = (1.0 + zeq)/(1.0 + zd); //matches wayne's paper, but not the code (1+zeq -> zeq)
  //fprintf(stderr,"paper eh 98 2\n");
#else
  y = (zeq)/(1.0 + zd);
#endif
  Gy = y*(-6.0*sqrt(1.0 + y) + (2.0 + 3.0*y)*log((sqrt(1.0 + y) + 1.0)/(sqrt(1.0 + y) - 1.0)));

  //eqn 14
  ab = 2.07*keq*s*pow(1.0 + Rd,-3.0/4.0)*Gy;

  //----------------------------------
  // Get CDM part of transfer function
  //----------------------------------

  //eqn 18
  f = 1.0/(1.0 + (k*s/5.4)*(k*s/5.4)*(k*s/5.4)*(k*s/5.4));

  //eqn 20
  C = 14.2/ac + 386.0/(1.0 + 69.9*pow(q,1.08));

  //eqn 19
  T0t = log(M_E + 1.8*bc*q)/(log(M_E + 1.8*bc*q) + C*q*q);

  //eqn 17
  C1bc = 14.2 + 386.0/(1.0 + 69.9*pow(q,1.08));
  T0t1bc = log(M_E + 1.8*bc*q)/(log(M_E + 1.8*bc*q) + C1bc*q*q);
  Tc = f*T0t1bc + (1.0 - f)*T0t;

  //-------------------------------------
  // Get baryon part of transfer function
  //-------------------------------------

  //eqn 24
  bb = 0.5 + omb/om0 + (3.0 - 2.0*omb/om0)*sqrt((17.2*om0*h*h)*(17.2*om0*h*h) + 1.0);

  //eqn 23
  bnode = 8.41*pow(om0*h*h,0.435);

  //eqn 22
  st = s/pow(1.0 + (bnode/k/s)*(bnode/k/s)*(bnode/k/s),1.0/3.0);

  //eqn 21
  C11 = 14.2 + 386.0/(1.0 + 69.9*pow(q,1.08));
  T0t11 = log(M_E + 1.8*q)/(log(M_E + 1.8*q) + C11*q*q);
  Tb = (T0t11/(1.0 + (k*s/5.2)*(k*s/5.2))
	+ ab/(1.0 + (bb/k/s)*(bb/k/s)*(bb/k/s))/exp(pow(k/ksilk,1.4)))*sin(k*st)/(k*st);

  //------------------------
  // total transfer function
  //------------------------
  Tk = omb/om0*Tb + omc/om0*Tc;

  return Tk;
}

double transfunct_eh98_smooth(double kin)
{
  //vars
  double k,Tk;
  double omb,om0,omc,h;
  double theta2p7,s,q;
  double Gamma,alphaGamma,L0,C0;

  //get cosmoparms
  omb = cosmoData.OmegaB;
  om0 = cosmoData.OmegaM;
  omc = cosmoData.OmegaM - cosmoData.OmegaB;
  h = cosmoData.h;

  //convert k from hMpc^-1 to Mpc^-1
  k = kin*h;

  //-----------
  //input parms
  //-----------
  theta2p7 = TCMB/2.7;

  //eqn 26
  s = 44.5*log(9.83/om0/h/h)/sqrt(1.0 + 10.0*pow(omb*h*h,0.75));

  //eqn 31
  alphaGamma = 1.0 - 0.328*log(431.0*om0*h*h)*omb/om0 + 0.38*log(22.3*om0*h*h)*(omb/om0)*(omb/om0);

  //eqn 30
  Gamma = om0*h*(alphaGamma + (1.0 - alphaGamma)/(1.0 + pow(0.43*k*s,4.0)));

  //eqn 28
  q = kin*theta2p7*theta2p7/Gamma;

  //eqns 29
  C0 = 14.2 + 731.0/(1.0 + 62.5*q);
  L0 = log(2.0*exp(1.0) + 1.8*q);
  Tk = L0/(L0 + C0*q*q);

  return Tk;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc