text: #define PJ_LIB__
  
  #include <errno.h>
  #include <math.h>
  
  #include "proj.h"
  #include "proj_internal.h"
  
  PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
  
  
  # define C1 .16666666666666666666
  # define C2 .00833333333333333333
  # define C3 .04166666666666666666
  # define C4 .33333333333333333333
  # define C5 .06666666666666666666
  
  
  namespace { // anonymous namespace
  struct cass_data {
      double *en;
      double m0;
      bool hyperbolic;
  };
  } // anonymous namespace
  
  
  
  static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
      PJ_XY xy = {0.0, 0.0};
      struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);
  
      const double sinphi = sin (lp.phi);
      const double cosphi = cos (lp.phi);
      const double M = pj_mlfn (lp.phi, sinphi, cosphi, Q->en);
  
      const double nu_square = 1./(1. - P->es * sinphi*sinphi);
      const double nu  = sqrt(nu_square);
      const double tanphi = tan(lp.phi);
      const double T = tanphi * tanphi;
      const double A = lp.lam * cosphi;
      const double C = P->es * (cosphi * cosphi) / (1 - P->es);
      const double A2 = A * A;
  
      xy.x = nu * A * (1. - A2 * T *
          (C1 - (8. - T + 8. * C) * A2 * C2));
      xy.y = M - Q->m0 + nu * tanphi * A2 *
          (.5 + (5. - T + 6. * C) * A2 * C3);
      if( Q->hyperbolic )
      {
          const double rho = nu_square * (1. - P->es) * nu;
          xy.y -= xy.y * xy.y * xy.y / (6 * rho * nu);
      }
  
      return xy;
  }
  
  
  static PJ_XY cass_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
      PJ_XY xy = {0.0, 0.0};
      xy.x  =  asin (cos (lp.phi) * sin (lp.lam));
      xy.y  =  atan2 (tan (lp.phi), cos (lp.lam)) - P->phi0;
      return xy;
  }
  
  
  static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
      PJ_LP lp = {0.0, 0.0};
      struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);
  
      const double phi1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
      const double tanphi1 = tan (phi1);
      const double T1 = tanphi1*tanphi1;
      const double sinphi1 = sin (phi1);
      const double nu1_square = 1. / (1. - P->es * sinphi1 * sinphi1);
      const double nu1 = sqrt (nu1_square);
      const double rho1 = nu1_square * (1. - P->es) * nu1;
      const double D  = xy.x / nu1;
      const double D2 = D * D;
      lp.phi = phi1 - (nu1 * tanphi1 / rho1) * D2 *
          (.5 - (1. + 3. * T1) * D2 * C3);
      lp.lam = D * (1. + T1 * D2 *
          (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);
  
      if( Q->hyperbolic )
      {
          // EPSG guidance note 7-2 suggests a custom approximation for the
          // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
          // generic inversion method
          lp = pj_generic_inverse_2d(xy, P, lp);
      }
  
      return lp;
  }
  
  
  static PJ_LP cass_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
      PJ_LP lp = {0.0,0.0};
      double dd;
      lp.phi = asin(sin(dd = xy.y + P->phi0) * cos(xy.x));
      lp.lam = atan2(tan(xy.x), cos(dd));
      return lp;
  }
  
  static PJ *destructor (PJ *P, int errlev) {                        /* Destructor */
      if (nullptr==P)
          return nullptr;
  
      if (nullptr==P->opaque)
          return pj_default_destructor (P, errlev);
  
      free (static_cast<struct cass_data*>(P->opaque)->en);
      return pj_default_destructor (P, errlev);
  }
  
  
  
  PJ *PROJECTION(cass) {
  
      /* Spheroidal? */
      if (0==P->es) {
          P->inv = cass_s_inverse;
          P->fwd = cass_s_forward;
          return P;
      }
  
      /* otherwise it's ellipsoidal */
      auto Q = static_cast<struct cass_data*>(calloc (1, sizeof (struct cass_data)));
      P->opaque = Q;
      if (nullptr==P->opaque)
          return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
      P->destructor = destructor;
  
      Q->en = pj_enfn (P->es);
      if (nullptr==Q->en)
          return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
  
      Q->m0 = pj_mlfn (P->phi0,  sin (P->phi0),  cos (P->phi0), Q->en);
      if (pj_param_exists(P->params, "hyperbolic"))
          Q->hyperbolic = true;
      P->inv = cass_e_inverse;
      P->fwd = cass_e_forward;
  
      return P;
  }
  