Civil twilight computation is using raw solar angle; shouldn't

[ericberman]

Needs to account for refraction, etc. Reverse engineer https://gml.noaa.gov/grad/solcalc/azel.html?

[ericberman]

Looks like this code does the trick of reverse engineering the spreadsheet at https://gml.noaa.gov/grad/solcalc/NOAA_Solar_Calculations_day.xls. Need to port to iOS and Android

        static public double calcSolarAngle(double lat, double lon, double JD, double minutes)
        {
            double julianCentury = calcTimeJulianCent(JD + minutes / 1440.0);

            double sunDeclinationRad = degToRad(calcSunDeclination(julianCentury));
            double latRad = degToRad(lat);

            double eqOfTime = calcEquationOfTime(julianCentury);

            double trueSolarTimeMin =  (minutes + eqOfTime + 4 *lon) % 1440;
            double hourAngleDeg = trueSolarTimeMin / 4 < 0 ? trueSolarTimeMin / 4 + 180 : trueSolarTimeMin / 4 - 180;
            double zenith = radToDeg(Math.Acos(Math.Sin(latRad) * Math.Sin(sunDeclinationRad) + Math.Cos(latRad) * Math.Cos(sunDeclinationRad) * Math.Cos(degToRad(hourAngleDeg))));
            double solarElevation = 90 - zenith;
            double atmRefractionDeg = solarElevation > 85 ? 0 :
                (solarElevation > 5 ? 58.1 / Math.Tan(degToRad(solarElevation)) - 0.07 / Math.Pow(Math.Tan(degToRad(solarElevation)), 3) + 0.000086 / Math.Pow(Math.Tan(degToRad(solarElevation)), 5) :
                solarElevation > -0.575 ? 1735 + solarElevation * (-518.2 + solarElevation * (103.4 + solarElevation * (-12.79 + solarElevation * 0.711))) : -20.772 / Math.Tan(degToRad(solarElevation))) / 3600;
            return solarElevation + atmRefractionDeg;
        }