ProJ in VEX

#include <math.h>
float theta, phi;
float r = chf("radius");

if(chi("reverse")) {
    theta = asin(@P.z/r);
    phi = atan2(@P.y,@P.x);
    float lat = degrees(theta);
    float lon = degrees(phi);
    @P.x = lon;
    @P.z = lat*-1;
    @P.y = 0;
} else {

    
    float lat1 = vector(getbbox_max(0)).z;
    float lat2 = vector(getbbox_min(0)).z;
    float lon1 = vector(getbbox_min(0)).x;
    float lon2 = vector(getbbox_max(0)).x;
    
    
    theta = PI/2 - radians(@P.z*-1);
    phi = radians(@P.x);
    
    float dLat = lat2 * PI / 180 - lat1 * PI / 180;
    float dLon = lon1 * PI / 180 - lon1 * PI / 180;
    
    float a = sin(dLat/2) * sin(dLat/2) + cos(lat1 
    * PI / 180) * cos(lat2 * PI / 180) 
    * sin(dLon/2) * sin(dLon/2);
    float c = 2 * atan2(sqrt(a), sqrt(1-a))*r;
    @d = r*c;
    float x = (r+@P.y*c) * sin(theta) * cos(phi);
    float y = (r+@P.y*c) * sin(theta) * sin(phi);
    float z = (r+@P.y*c) * cos(theta);
    v@P = set(x,y,z);
}

//def gps_to_ecef_custom(lat, lon, alt):
//WGS84 to ECEF

#include <math.h>
float rad_lat = (@P.z) * (PI / 180.0);
float rad_lon = @P.x * (PI / 180.0);

float a = 6378.1370;
float finv = 298.257223563;
float f = 1 / finv;

float e2 = 1 - (1 - f) * (1 - f);
float v = a / sqrt(1 - e2 * sin(rad_lat) * sin(rad_lat));

@P.x = (v + @P.y) * cos(rad_lat) * cos(rad_lon);
@P.z = (v + @P.y) * cos(rad_lat) * sin(rad_lon) *-1;
@P.y = (v * (1 - e2) + @P.y) * sin(rad_lat) *-1;

//Equal Earth Projection
//https://en.wikipedia.org/wiki/Equal_Earth_projection

float A1 = 1.340264;
float A2 = -0.081106;
float A3 = 0.000893;
float A4 = 0.003796;
float M = sqrt(3) / 2.0;
float paramLatSq, paramLatPow6, paramLat;

paramLat = asin(M * sin(radians(@P.z)));
paramLatSq = paramLat * paramLat;
paramLatPow6 = paramLatSq * paramLatSq * paramLatSq;

@P.x = radians(@P.x) * cos(paramLat)
        / (M * (A1 + 3 * A2 * paramLatSq + paramLatPow6 * (7 * A3 + 9 * A4 * paramLatSq)));
@P.z = paramLat * (A1 + A2 * paramLatSq + paramLatPow6 * (A3 + A4 * paramLatSq));

//Gall Peters 
//https://en.wikipedia.org/wiki/Gall%E2%80%93Peters_projection


#include <math.h>
float R = 6378.1370;
vector pos = @P;

@P.x = R*PI*@P.x*cos(radians(45))/180;
@P.z = R*sin(radians(@P.z)) / cos(radians(45));

//@P.x = R*PI*pos.x / 180 * sqrt(2);
//@P.z = R * sqrt(2) * sin(radians(pos.z));

//@P.x = R*pos.x;
//@P.y = 2*R*sin(pos.z);

//https://desktop.arcgis.com/de/arcmap/latest/map/projections/craster-parabolic.htm

float XM = 0.97720502380583984317;
float RXM = 1.02332670794648848847;
float YM = 3.06998012383946546542;
float RYM = 0.32573500793527994772;
float THIRD = 0.333333333333333333;

float lpphi  = radians(@P.z);
float lplam = radians(@P.x);
lpphi *= THIRD;
@P.x = XM * lplam * (2. * cos((lpphi + lpphi)) - 1.);
@P.z = YM * sin(lpphi);
@P *= chf("radius");

//https://www.youtube.com/watch?v=D3tdW9l1690
//https://en.wikipedia.org/wiki/Euler_spiral

//INIT START
float projectionLatitude = radians(ch("eulerLat0"));
float projectionLatitude1 = radians(ch("eulerLat1"));
float projectionLatitude2 = radians(ch("eulerLat2"));

@del = 0.5 * (projectionLatitude2 - projectionLatitude1);
float sig = 0.5 * (projectionLatitude2 + projectionLatitude1);

@n = sin(sig) * sin((@del)) / @del;
@del *= 0.5;
@rho_c = @del / (tan(@del) * tan(sig)) + sig;
@rho_0 = @rho_c - projectionLatitude;
//INIT END

float lplam = radians(@P.x);
float lpphi = radians(@P.z*-1);

float rho = @rho_c - lpphi;
@P.x = rho * sin(lplam*=@n+radians(chf("radius")));
@P.z = @rho_0 - rho * cos(lplam);

float XM = 0.97720502380583984317;
float RXM = 1.02332670794648848847;
float YM = 3.06998012383946546542;
float RYM = 0.32573500793527994772;
float THIRD = 0.333333333333333333;
float lpphi = radians(@P.z*-1);
float lplam = radians(@P.x);
lpphi *= THIRD;
@P.x = XM * lplam * (2.0 * cos(lpphi + lpphi) - 1.0)*1.0;
@P.z = YM * sin(lpphi)*1.0;

#include <math.h>

float scaleFactor = chf("scaleFactor");
float phi = radians(@P.z*-1);
float lam = radians(@P.x);
float maxLatitude = radians(85);
float minLatitude = radians(-85);
float e = 2.71828;

float tsfn(float phi; float sinphi; float e;) 
{
    sinphi *= e;
    return (tan(.5 * ((PI/2) - phi))
            / pow((1. - sinphi) / (1. + sinphi), .5 * e));
}
float __sinphi__ = sin(radians(phi));

if (chi("spherical")) {
    @P.x = scaleFactor * lam;
    if (phi > maxLatitude) {
        phi = maxLatitude;
    } else if (phi < minLatitude) {
        phi = minLatitude;
    }
    @P.z = scaleFactor * log(tan((PI/4) + 0.5 * phi));
} else {
    @P.x = scaleFactor * lam;
    @P.z = -scaleFactor * log( tsfn(phi, __sinphi__, e) );
}

#include <math.h>
//Eckert V https://proj.org/operations/projections/eck5.html
float XF = 0.44101277172455148219;
float RXF = 2.26750802723822639137;
float YF = 0.88202554344910296438;
float RYF = 1.13375401361911319568;

float lplam = radians(@P.x);
float lpphi = (radians(@P.z*-1));

@P.x = XF * (1. + cos(lpphi)) * lplam;
@P.z = YF * lpphi;

#include <math.h>

float w = 0.25;
float m = 1;
float rm;

if ((m = abs(m)) <= 0.) {
    printf("ERROR");
} else {
    m = 1.;
}
rm = 1.0 / m;
m /= w;

float cosphi, d;
float lpphi = radians(@P.z*-1);
float lplam = radians(@P.x);
d = sqrt(2. / (1. + (cosphi = cos(lpphi)) * cos(lplam *= w)));
//cosphi = cos(lpphi) * cos(lplam *= w);
@P.x = m * d * cosphi * sin(lplam) ;
@P.z = rm * d * sin(lpphi);

//https://en.wikipedia.org/wiki/Winkel_tripel_projection
float lplam = radians(@P.x);
float lpphi = radians(@P.z);

float c = 0.5 * lplam;
float d = acos(cos(lpphi) * cos(c));

if (d != 0) {
    @P.x = 2. * d * cos(lpphi) * sin(c) * (@P.z = 1. / sin(d));
    @P.z *= d * sin(lpphi);
} else {
    @P.x = @P.z = 0.0;
}
@P.x = (@P.x + lplam * 0.636619772367581343) * 0.5;
@P.z = (@P.z + lpphi) * 0.5;

#include <math.h>


float phi1 = acos(2./PI);


float cosphi1 = 2./PI;
int MAX_ITER = 10;
float LOOP_TOL = 1e-7;
float TWO_D_PI = 0.636619772367581343;

float k, V;
int i;


float lpphi = radians(@P.z);
float lplam = radians(@P.x);
@P.z = lpphi * TWO_D_PI;
k = PI * sin(lpphi);
lpphi *= 1.8;

for (i = MAX_ITER; i > 0; --i) {
    lpphi -= V = (lpphi + sin(lpphi) - k) / (1. + cos(lpphi));
    if (abs(V) < LOOP_TOL)
        break;
}

@P.x = 0.5 * lplam * (cos(lpphi) + cosphi1);
@P.z = (PI/4) * (sin(lpphi) + @P.z);

float phi1 = radians(50. + 28. / 60.);
float cosphi1 = cos(phi1);
float lplam = radians(@P.x);
float lpphi = radians(@P.z);

@P.x = .5 * lplam * (cosphi1 + cos(lpphi));
@P.z = lpphi;

float theta, ct, D;
float lpphi = radians(@P.z);
float lplam = radians(@P.x);

theta = asin(@P.z = 0.90630778703664996 * sin(lpphi));
@P.x = 2.66723 * (ct = cos(theta)) * sin(lplam /= 3.);
@P.z *= 1.24104 * (D = 1 / (sqrt(0.5 * (1 + ct * cos(lplam)))));
@P.x *= D;

//VanDerGrintenProjection
//https://en.wikipedia.org/wiki/Van_der_Grinten_projection
#include <math.h>

float TOL = 1.e-10;
float THIRD = .33333333333333333333;
float TWO_THRD = .66666666666666666666;
float C2_27 = .07407407407407407407;
float PI4_3 = 4.18879020478639098458;
float PISQ = 9.86960440108935861869;
float TPISQ = 19.73920880217871723738;
float HPISQ = 4.93480220054467930934;

float lplam = radians(@P.x);
float lpphi = radians(@P.z);

float  al, al2, g, g2, p2;
p2 = abs(lpphi / (PI/2));

if ((p2 - TOL) > 1.) printf("ERROR #F");
if (p2 > 1.)
        p2 = 1.;
if (abs(lpphi) <= TOL) {
        @P.x = lplam;
        @P.z = 0.;
} else if (abs(lplam) <= TOL || abs(p2 - 1.) < TOL) {
        @P.x = 0.;
        @P.z = PI * tan(.5 * asin(p2));
        if (lpphi < 0.) @P.z = -@P.z;
} else {
        al = .5 * abs(PI / lplam - lplam / PI);
        al2 = al * al;
        g = sqrt(1. - p2 * p2);
        g = g / (p2 + g - 1.);
        g2 = g * g;
        p2 = g * (2. / p2 - 1.);
        p2 = p2 * p2;
        @P.x = g - p2; g = p2 + al2;
        @P.x = PI * (al * @P.x + sqrt(al2 * @P.x * @P.x - g * (g2 - p2))) / g;
        if (lplam < 0.) @P.x = -@P.x;
        @P.z = abs(@P.x / PI);
        @P.z = 1. - @P.z * (@P.z + 2. * al);
        if (@P.z < -TOL) printf("ERROR #F");
        if (@P.z < 0.)
                @P.z = 0.;
        else
                @P.z = sqrt(@P.z) * (lpphi < 0. ? -PI : PI);
}

//URMFPSProjection
// * With the default parameter n, this is identical to Wagner I.

float C_x = 0.8773826753;
float C_y = 1.139753528477;
float n = 0.8660254037844386467637231707;// wag1
//float C_y = chf("radius");

float lplam = radians(@P.x);
float lpphi = radians(@P.z);

lpphi = asin(n * sin(lpphi));
@P.x = C_x * lplam * cos(lpphi);
@P.z = C_y * lpphi;

float projectionLatitude = radians(45);
float projectionLatitude1 = radians(chf("lat_1_35"));
float projectionLatitude2 = radians(chf("lat_2_60"));
float EPS = 1e-10;
float rho_c, rho_0, n;

// FIXME this will throw an exception with lat1 = lat2, which is the default.
float del = 0.5 * (projectionLatitude2 - projectionLatitude1);
float sig = 0.5 * (projectionLatitude2 + projectionLatitude1);

if (abs(del) < EPS || abs(sig) < EPS) {
    //throw new ProjectionException("-42");
}
n = sin(sig);
float cs = cos(del);
rho_c = n / cs + cs / n;
rho_0 = sqrt((rho_c - 2 * sin(projectionLatitude)) / n);

float lpphi = radians(@P.z*-1);
float lplam = radians(@P.x);

float rho = rho_c - lpphi;
@P.x = rho * sin(lplam *= n);
@P.z = rho_0 - rho * cos(lplam)*-1;