How to calculate sunrise and sunset times?

Here is a complete routine to that calculates the sunrise (optionally the sunset, see the code) in C++. It only requires latitude and longitude as input, and is accurate to within seconds for the civil sunrise time. This code is running in a UWP C++ app created with Visual Studio 2017. It includes a subroutine to fill the output minutes and seconds with zeros in the case of single digit results.

String^ toplatformstring(bool fill, long ll) {
    // convert ll to platform string
    Platform::String^ p_string;
    std::string doit;

    if (fill == false) {
        doit = std::to_string(ll); // convert, don't fill with zeros
    }
    else {
        //convert ll to std doit and fill with zeros
        std::stringstream ss;
        ss << std::setw(2) << std::setfill('0') << ll;
        doit = ss.str();
    }

    //convert doit to platform string
    char const *pchar = doit.c_str();
    std::string s_str = std::string(pchar);
    std::wstring wid_str = std::wstring(s_str.begin(), s_str.end());
    const wchar_t* w_char = wid_str.c_str();
    p_string = ref new Platform::String(w_char);

    return p_string;
}

//double la = 39.299236;  // baltimore
//double lo = -76.609383;
//double la = 37.0;  // SF California
//double lo = -122.0;
Platform::String^ MainPage::sunrise(double la, double lo) {


    /*double la = 39.300213;
    double lo = -76.610516;*/
    Platform::String^ sunrisetime;

    //// get year, month, day integers
    time_t rawtime;
    struct tm timeinfo;  // get date and time info
    time(&rawtime);
    localtime_s(&timeinfo, &rawtime);

    double xday = timeinfo.tm_mday;
    double xmonth = timeinfo.tm_mon;
    xmonth = xmonth + 1; // correct for origin 0
    //textblockc->Text = xmonth.ToString();  // for debugging
    double xyear = timeinfo.tm_year;
    //double dayofyear = timeinfo.tm_yday; // day of year also

    // calculate the day of the year
    //  N1 = floor(275 * month / 9)
    double xxN1 = floor(275 * xmonth / 9);
    //  N2 = floor((month + 9) / 12)
    double xxN2 = floor((xmonth + 9) / 12);
    //  N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
    double xxN3 = (1 + floor((xyear - 4 * floor(xyear / 4) + 2) / 3));
    //  N = N1 - (N2 * N3) + day - 30
    double day = xxN1 - (xxN2 * xxN3) + xday - 30;

    double zenith = 90.83333333333333;
    double D2R = M_PI / 180;
    double R2D = 180 / M_PI;

    // convert the longitude to hour value and calculate an approximate time
    double lnHour = lo / 15;
    double t;
    //if (sunrise) {
    t = day + ((6 - lnHour) / 24);
    //} else {
    //t = day + ((18 - lnHour) / 24);
    //};

    //calculate the Sun's mean anomaly
    double M = (0.9856 * t) - 3.289;

    //calculate the Sun's true longitude
    double L = M + (1.916 * sin(M * D2R)) + (0.020 * sin(2 * M * D2R)) + 282.634;
    if (L > 360) {
        L = L - 360;
    }
    else if (L < 0) {
        L = L + 360;
    };

    //calculate the Sun's right ascension
    double RA = R2D * atan(0.91764 * tan(L * D2R));
    if (RA > 360) {
        RA = RA - 360;
    }
    else if (RA < 0) {
        RA = RA + 360;
    };

    //right ascension value needs to be in the same qua
    double Lquadrant = (floor(L / (90))) * 90;
    double RAquadrant = (floor(RA / 90)) * 90;
    RA = RA + (Lquadrant - RAquadrant);

    //right ascension value needs to be converted into hours
    RA = RA / 15;

    //calculate the Sun's declination
    double sinDec = 0.39782 * sin(L * D2R);
    double cosDec = cos(asin(sinDec));

    //calculate the Sun's local hour angle
    double cosH = (cos(zenith * D2R) - (sinDec * sin(la * D2R))) / (cosDec * cos(la * D2R));
    double H;
    //if (sunrise) {
    H = 360 - R2D * acos(cosH);
    //} else {
    //H = R2D * Math.acos(cosH)
    //};
    H = H / 15;

    //calculate local mean time of rising/setting
    double T = H + RA - (0.06571 * t) - 6.622;

    //adjust back to UTC
    double UT = T - lnHour;
    if (UT > 24) {
        UT = UT - 24;
    }
    else if (UT < 0) {
        UT = UT + 24;
    }

    //convert UT value to local time zone of latitude/longitude
    int offset = (int)(lo / 15); // estimate utc correction
    double localT = UT + offset; // -5 for baltimore

                                 //convert to seconds
    int seconds = (int)(localT * 3600);

    long sec = seconds % 60;
    long minutes = seconds % 3600 / 60;
    long hours = seconds % 86400 / 3600;
    hours = hours % 12;

    Platform::String^ ssec = toplatformstring(true, sec);
    Platform::String^ mminutes = toplatformstring(true, minutes);
    Platform::String^ hhours = toplatformstring(false, hours);


    sunrisetime = hhours + ":" + mminutes + ":" + ssec;
    return sunrisetime;
}