/***********************************************/
/*           Copyright (c) 2025 Belmu          */
/*             All Rights Reserved             */
/***********************************************/

/*
    [References]:
        Polyanskiy, M. (2017). Ciddor 1996 - air. https://github.com/polyanskiy/refractiveindex.info-scripts/blob/master/scripts/Ciddor%201996%20-%20air.py
        SCHOTT. (2021). Downloads for Optical Glass. https://www.schott.com/en-be/products/optical-glass-p1000267/downloads
*/

struct SellmeierCoefficients {
    vec3 B, C;
};

const SellmeierCoefficients N_AIR = SellmeierCoefficients(vec3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 0.0));

const SellmeierCoefficients N_BK7HT = SellmeierCoefficients(vec3(1.039612120, 0.231792344, 1.010469450), vec3(0.00600069867, 0.0200179140, 103.56065300));
const SellmeierCoefficients N_LAF2  = SellmeierCoefficients(vec3(1.809842270, 0.157295550, 1.093003700), vec3(0.01017116200, 0.0442431765, 100.68774800));
const SellmeierCoefficients N_SF8   = SellmeierCoefficients(vec3(1.550758120, 0.209816918, 1.462054910), vec3(0.01143383400, 0.0582725652, 133.24165000));
const SellmeierCoefficients F2HT    = SellmeierCoefficients(vec3(1.345333590, 0.209073176, 0.937357162), vec3(0.00997743900, 0.0470450767, 111.88676400));

const SellmeierCoefficients N_BASF2  = SellmeierCoefficients(vec3(1.536520810, 0.156971102, 1.301968150), vec3(0.010843573, 0.0562278762, 131.3397000));
const SellmeierCoefficients N_BAF10  = SellmeierCoefficients(vec3(1.585149500, 0.143559385, 1.085212690), vec3(0.009266813, 0.0424489805, 105.6135730));
const SellmeierCoefficients N_BAF4   = SellmeierCoefficients(vec3(1.420563280, 0.102721269, 1.143809760), vec3(0.009420154, 0.0531087291, 110.2788560));
const SellmeierCoefficients N_SSK5   = SellmeierCoefficients(vec3(1.592226590, 0.103520774, 1.051740160), vec3(0.009202846, 0.0423530072, 106.9273740));
const SellmeierCoefficients N_K5     = SellmeierCoefficients(vec3(1.085118330, 0.199562005, 0.930511663), vec3(0.006610995, 0.0241108660, 111.9827770));
const SellmeierCoefficients N_SF5    = SellmeierCoefficients(vec3(1.524818890, 0.187085527, 1.427290150), vec3(0.011254756, 0.0588995392, 129.1416750));
const SellmeierCoefficients N_SF64   = SellmeierCoefficients(vec3(1.591637620, 0.219908428, 1.469293150), vec3(0.011862343, 0.0594585499, 133.3107620));
const SellmeierCoefficients SF1      = SellmeierCoefficients(vec3(1.559129230, 0.284246288, 0.968842926), vec3(0.012148100, 0.0534549042, 112.1748090));
const SellmeierCoefficients SF2      = SellmeierCoefficients(vec3(1.403018210, 0.231767504, 0.939056586), vec3(0.010579547, 0.0493226978, 112.4059550));
const SellmeierCoefficients SF6      = SellmeierCoefficients(vec3(1.724484820, 0.390104889, 1.045728580), vec3(0.013487195, 0.0569318095, 118.5571850));
const SellmeierCoefficients F5       = SellmeierCoefficients(vec3(1.310446300, 0.196034260, 0.966129770), vec3(0.009586330, 0.0457627627, 115.0118830));
const SellmeierCoefficients LF7      = SellmeierCoefficients(vec3(1.268797310, 0.156697051, 0.887676301), vec3(0.009195539, 0.0447256764, 106.8203610));
const SellmeierCoefficients N_F2     = SellmeierCoefficients(vec3(1.397570370, 0.159201403, 1.268654300), vec3(0.009959061, 0.0546931752, 119.2483460));
const SellmeierCoefficients N_LAK8   = SellmeierCoefficients(vec3(1.331831670, 0.546623206, 1.190840150), vec3(0.006200239, 0.0216465439, 082.5827736));
const SellmeierCoefficients N_LAK9   = SellmeierCoefficients(vec3(1.462319050, 0.344399589, 1.155083720), vec3(0.007242702, 0.0243353131, 085.4686868));
const SellmeierCoefficients N_LAK21  = SellmeierCoefficients(vec3(1.227181160, 0.420783743, 1.012848430), vec3(0.006020757, 0.0196862889, 088.4370099));
const SellmeierCoefficients N_LAK33B = SellmeierCoefficients(vec3(1.422886010, 0.593661336, 1.161352600), vec3(0.006702835, 0.0219416210, 080.7407701));
const SellmeierCoefficients N_KZFS11 = SellmeierCoefficients(vec3(1.332224500, 0.289241610, 1.151617340), vec3(0.008402985, 0.0344239720, 088.4310532));
const SellmeierCoefficients P_SK57   = SellmeierCoefficients(vec3(1.310534140, 0.169376189, 1.109877140), vec3(0.007408772, 0.0254563489, 107.7510870));

const SellmeierCoefficients S_BAH11  = SellmeierCoefficients(vec3(1.571388600, 0.147869313, 1.280928460), vec3(0.009108079, 0.0402401684, 130.3993670));
const SellmeierCoefficients S_LAH65V = SellmeierCoefficients(vec3(1.814190340, 0.361376301, 1.327294840), vec3(0.008749350, 0.0318352836, 091.3406898));
const SellmeierCoefficients S_LAM3   = SellmeierCoefficients(vec3(1.642587130, 0.239634610, 1.224830260), vec3(0.008682460, 0.0351226242, 116.6043690));

float sellmeier(SellmeierCoefficients coefficients, float lambda) {
    if (coefficients == N_AIR) return airIOR;

    lambda *= 1e-3; // To micrometers
    float lambdaSq = lambda * lambda;

    vec3 eta = (coefficients.B * lambdaSq) / (vec3(lambdaSq) - coefficients.C);
    return sqrt(1.0 + eta.x + eta.y + eta.z);
}

float waterDispersion(float lambda) {
    lambda *= 1e-3; // To micrometers
    float lambdaSq = lambda * lambda;
    
    float dispersion  = 1.0;
          dispersion += ((5.689093832e-1 * lambdaSq) / (lambdaSq - 5.110301794e-3));
          dispersion += ((1.719708856e-1 * lambdaSq) / (lambdaSq - 1.825180155e-2));
          dispersion += ((2.062501582e-2 * lambdaSq) / (lambdaSq - 2.624158904e-2));
          dispersion += ((1.123965424e-1 * lambdaSq) / (lambdaSq - 1.067505178e1 ));

    return dispersion;
}

const float absoluteZero = 273.15;

float compressibility(float T, float p, float xw) {
    T -= absoluteZero;

    const float a0 =  1.58123e-6;
    const float a1 = -2.9331e-8;
    const float a2 =  1.1043e-10;
    const float b0 =  5.707e-6;
    const float b1 = -2.051e-8;
    const float c0 =  1.9898e-4;
    const float c1 = -2.376e-6;
    const float d  =  1.83e-11;
    const float e  = -0.765e-8;
    return 1.0 - (p / T) * (a0 + a1 * T + a2 * T * T + (b0 + b1 * T) * xw + (c0 + c1 * T) * xw * xw) + pow2(p / T) * (d + e * xw * xw);
}

float airDispersion(float lambda, float T, float p, float h, float xc) {
    lambda *= 1e-3; // To micrometers

    T += absoluteZero;

    const float R     =  8.314510;
    const float k0    =  238.0185;
    const float k1    =  5792105.0;
    const float k2    =  57.362;
    const float k3    =  167917.0;
    const float w0    =  295.235;
    const float w1    =  2.6422;
    const float w2    = -0.032380;
    const float w3    =  0.004028;
    const float A     =  1.2378847e-5;
    const float B     = -1.9121316e-2;
    const float C     =  33.93711047;
    const float D     =  6.3431645e3;
    const float alpha =  1.00062;
    const float beta  =  3.14e-8;
    const float gamma =  5.6e-7;

    float sigma  = 1.0 / lambda;
    float sigma2 = sigma * sigma;

    float svp = 0.0;

    // Saturation vapor pressure of water vapor in air at temperature T
    if (T >= 0.0) {
        svp = exp(A * pow2(T) + B * T + C + D / T);
    } else {
        svp = pow(10.0, -2663.5 / T + 12.537);
    }
    
    float f = alpha + beta * p + gamma * pow2(T); // Enhancement factor of water vapor in air
    
    float xw = f * h * svp / p; // Molar fraction of water vapor in moist air

    // Refractive index of standard air at 15 °C, 101325 Pa, 0% humidity, 450 ppm CO2
    float nas = 1.0 + (k1 / (k0 - sigma2) + k3 / (k2 - sigma2)) * 1e-8;
    // Refractive index of standard air at 15 °C, 101325 Pa, 0% humidity, xc ppm CO2
    float naxs = 1.0 + (nas - 1.0) * (1.0 + 0.534e-6 * (xc - 450.0));
    // Refractive index of water vapor at standard conditions (20 °C, 1333 Pa)
    float nws = 1.0 + 1.022 * (w0 + w1 * sigma2 + w2 * pow4(sigma) + w3 * pow6(sigma)) * 1e-8;

          float Ma = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400.0)); // Molar mass of dry air
    const float Mw = 0.018015;                                    // Molar mass of water vapor

    float Za = compressibility(288.15, 101325.0, 0.0); // Compressibility of dry air
    float Zw = compressibility(293.15, 1333.0  , 1.0); // Compressibility of pure water vapor

    float paxs = 101325.0 * Ma / (Za * R * 288.15); // Density of standard air
    float pws  = 1333.0   * Mw / (Zw * R * 293.15); // Density of standard water vapor

    float pa = p * Ma / (compressibility(T, p, xw) * R * T) * (1.0 - xw); // Density of the dry component of the moist air 
    float pw = p * Mw / (compressibility(T, p, xw) * R * T) * xw;         // Density of the water vapor component

    return 1.0 + (pa / paxs) * (naxs - 1.0) + (pw / pws) * (nws - 1.0);
}
