11

PHP の date_sun_infodate_sunrisedate_sunset 関数は、どのように計算しているのでしょうか。 関数によって算出された値の精度を知りたいので、ご存知でしたら教えて下さい。

1 件の回答 1

16

date_sunrisedate_sunsetdate.cphp_do_date_sunrise_sunset関数を使います。date.ctimelibライブラリを使っています。

static void php_do_date_sunrise_sunset(INTERNAL_FUNCTION_PARAMETERS, int calc_sunset)
{
    double latitude = 0.0, longitude = 0.0, zenith = 0.0, gmt_offset = 0, altitude;
    double h_rise, h_set, N;
    timelib_sll rise, set, transit;
    long time, retformat = 0;
    int             rs;
    timelib_time   *t;
    timelib_tzinfo *tzi;
    char           *retstr;

    if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "l|ldddd", &time, &retformat, &latitude, &longitude, &zenith, &gmt_offset) == FAILURE) {
        RETURN_FALSE;
    }

    switch (ZEND_NUM_ARGS()) {
        case 1:
            retformat = SUNFUNCS_RET_STRING;
        case 2:
            latitude = INI_FLT("date.default_latitude");
        case 3:
            longitude = INI_FLT("date.default_longitude");
        case 4:
            if (calc_sunset) {
                zenith = INI_FLT("date.sunset_zenith");
            } else {
                zenith = INI_FLT("date.sunrise_zenith");
            }
        case 5:
        case 6:
            break;
        default:
            php_error_docref(NULL TSRMLS_CC, E_WARNING, "invalid format");
            RETURN_FALSE;
            break;
    }
    if (retformat != SUNFUNCS_RET_TIMESTAMP &&
        retformat != SUNFUNCS_RET_STRING &&
        retformat != SUNFUNCS_RET_DOUBLE)
    {
        php_error_docref(NULL TSRMLS_CC, E_WARNING, "Wrong return format given, pick one of SUNFUNCS_RET_TIMESTAMP, SUNFUNCS_RET_STRING or SUNFUNCS_RET_DOUBLE");
        RETURN_FALSE;
    }
    altitude = 90 - zenith;

    /* Initialize time struct */
    t = timelib_time_ctor();
    tzi = get_timezone_info(TSRMLS_C);
    t->tz_info = tzi;
    t->zone_type = TIMELIB_ZONETYPE_ID;

    if (ZEND_NUM_ARGS() <= 5) {
        gmt_offset = timelib_get_current_offset(t) / 3600;
    }

    timelib_unixtime2local(t, time);
    rs = timelib_astro_rise_set_altitude(t, longitude, latitude, altitude, 1, &h_rise, &h_set, &rise, &set, &transit);
    timelib_time_dtor(t);

    if (rs != 0) {
        RETURN_FALSE;
    }

    if (retformat == SUNFUNCS_RET_TIMESTAMP) {
        RETURN_LONG(calc_sunset ? set : rise);
    }
    N = (calc_sunset ? h_set : h_rise) + gmt_offset;

    if (N > 24 || N < 0) {
        N -= floor(N / 24) * 24;
    }

    switch (retformat) {
        case SUNFUNCS_RET_STRING:
            spprintf(&retstr, 0, "%02d:%02d", (int) N, (int) (60 * (N - (int) N)));
            RETURN_STRINGL(retstr, 5, 0);
            break;
        case SUNFUNCS_RET_DOUBLE:
            RETURN_DOUBLE(N);
            break;
    }
}

timelibtimelib_astro_rise_set_altitude関数に英語のコメントがありますが、ソースを見ると下記の計算を利用しているらしいです:

日の出計算

int timelib_astro_rise_set_altitude(timelib_time *t_loc, double lon, double lat, double altit, int upper_limb, double *h_rise, double *h_set, timelib_sll *ts_rise, timelib_sll *ts_set, timelib_sll *ts_transit)
{
    double  d,  /* Days since 2000 Jan 0.0 (negative before) */
    sr,         /* Solar distance, astronomical units */
    sRA,        /* Sun's Right Ascension */
    sdec,       /* Sun's declination */
    sradius,    /* Sun's apparent radius */
    t,          /* Diurnal arc */
    tsouth,     /* Time when Sun is at south */
    sidtime;    /* Local sidereal time */
    timelib_time *t_utc;
    timelib_sll   timestamp, old_sse;

    int rc = 0; /* Return cde from function - usually 0 */

    /* Normalize time */
    old_sse = t_loc->sse;
    t_loc->h = 12;
    t_loc->i = t_loc->s = 0;
    timelib_update_ts(t_loc, NULL);

    /* Calculate TS belonging to UTC 00:00 of the current day */
    t_utc = timelib_time_ctor();
    t_utc->y = t_loc->y;
    t_utc->m = t_loc->m;
    t_utc->d = t_loc->d;
    t_utc->h = t_utc->i = t_utc->s = 0;
    timelib_update_ts(t_utc, NULL);

    /* Compute d of 12h local mean solar time */
    timestamp = t_loc->sse;
    d = timelib_ts_to_juliandate(timestamp) - lon/360.0;

    /* Compute local sidereal time of this moment */
    sidtime = astro_revolution(astro_GMST0(d) + 180.0 + lon);

    /* Compute Sun's RA + Decl at this moment */
    astro_sun_RA_dec( d, &sRA, &sdec, &sr );

    /* Compute time when Sun is at south - in hours UT */
    tsouth = 12.0 - astro_rev180(sidtime - sRA) / 15.0;

    /* Compute the Sun's apparent radius, degrees */
    sradius = 0.2666 / sr;

    /* Do correction to upper limb, if necessary */
    if (upper_limb) {
        altit -= sradius;
    }

    /* Compute the diurnal arc that the Sun traverses to reach */
    /* the specified altitude altit: */
    {
        double cost;
        cost = (sind(altit) - sind(lat) * sind(sdec)) / (cosd(lat) * cosd(sdec));
        *ts_transit = t_utc->sse + (tsouth * 3600);
        if (cost >= 1.0) {
            rc = -1;
            t = 0.0;       /* Sun always below altit */

            *ts_rise = *ts_set = t_utc->sse + (tsouth * 3600);
        } else if (cost <= -1.0) {
            rc = +1;
            t = 12.0;      /* Sun always above altit */

            *ts_rise = t_loc->sse - (12 * 3600);
            *ts_set  = t_loc->sse + (12 * 3600);
        } else {
            t = acosd(cost) / 15.0;   /* The diurnal arc, hours */

            /* Store rise and set times - as Unix Timestamp */
            *ts_rise = ((tsouth - t) * 3600) + t_utc->sse;
            *ts_set  = ((tsouth + t) * 3600) + t_utc->sse;

            *h_rise = (tsouth - t);
            *h_set  = (tsouth + t);
        }
    }

    /* Kill temporary time and restore original sse */
    timelib_time_dtor(t_utc);
    t_loc->sse = old_sse;

    return rc;
}

date_sun_infotimelib関数を利用しています。

この回答はcmbuckleyStack Overflow回答からの情報です。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。