最近在玩 DuckDB，根据看过的番想到了一个数据集，于是用 SQL 搞了点花活
看看有没有佬友能发现 SQL 中二次元的部分
没有缩进是因为 Web Shell 对换行的 SQL 支持有点问题，复制过来懒得改了
LOAD icu;
SET VARIABLE min_oppositions = 2;
SET VARIABLE min_observation_years = 2;
SET VARIABLE mpcorb_url = 'https://www.minorplanetcenter.net/iau/MPCORB/MPCORB.DAT.gz';
SET VARIABLE observation_utc = date_trunc('second', current_timestamp AT TIME ZONE 'UTC');
SET VARIABLE encoded_observation_utc = replace( strftime( '%Y-%m-%dT%H:%M:%S', getvariable('observation_utc') ), ':', '%3A' );
SET VARIABLE miriade_url = 'https://ssp.imcce.fr/webservices/miriade/api/ephemcc.php?-name=p%3AEarth&-ep=' || getvariable('encoded_observation_utc') || '&-nbd=1&-tscale=UTC&-observer=%40sun&-theory=INPOP&-teph=4&-tcoor=2&-rplane=2&-mime=json&-output=--jd%2C--ttmutc&-from=duckdb';
CREATE TEMP TABLE earth_ephemeris AS WITH response AS ( SELECT d.Date::DOUBLE AS jd_ut, d."TT-UTC"::DOUBLE AS tt_minus_ut_seconds, d.px::DOUBLE AS earth_x_au, d.py::DOUBLE AS earth_y_au, d.pz::DOUBLE AS earth_z_au FROM read_json_auto(getvariable('miriade_url')) AS payload, UNNEST(payload.data) AS item(d) ), time_scales AS ( SELECT *, radians( 357.53 + 0.9856003 * (jd_ut - 2451545.0) ) AS earth_mean_anomaly_rad FROM response ) SELECT getvariable('observation_utc') AS observation_utc, jd_ut, tt_minus_ut_seconds + 0.001657 * sin(earth_mean_anomaly_rad) + 0.000022 * sin(2.0 * earth_mean_anomaly_rad) AS tdb_minus_ut_seconds, earth_x_au, earth_y_au, earth_z_au FROM time_scales;
CREATE TEMP MACRO packed_value(packed_character) AS ( CASE WHEN packed_character BETWEEN '0' AND '9' THEN ascii(packed_character) - ascii('0') ELSE ascii(packed_character) - ascii('A') + 10 END );
CREATE TEMP TABLE asteroid_distances AS WITH source_orbits AS ( SELECT trim(substr(line, 1, 7)) AS packed_designation, try_cast(trim(substr(line, 9, 5)) AS DOUBLE) AS absolute_magnitude_h, trim(substr(line, 21, 5)) AS packed_epoch, try_cast(trim(substr(line, 27, 9)) AS DOUBLE) AS mean_anomaly_deg, try_cast(trim(substr(line, 38, 9)) AS DOUBLE) AS argument_of_perihelion_deg, try_cast(trim(substr(line, 49, 9)) AS DOUBLE) AS ascending_node_deg, try_cast(trim(substr(line, 60, 9)) AS DOUBLE) AS inclination_deg, try_cast(trim(substr(line, 71, 9)) AS DOUBLE) AS eccentricity, try_cast(trim(substr(line, 81, 11)) AS DOUBLE) AS mean_daily_motion_deg, try_cast(trim(substr(line, 93, 11)) AS DOUBLE) AS semi_major_axis_au, try_cast(trim(substr(line, 118, 5)) AS INTEGER) AS observation_count, try_cast(trim(substr(line, 124, 3)) AS INTEGER) AS opposition_count, trim(substr(line, 128, 9)) AS observation_arc, trim(substr(line, 167, 28)) AS readable_designation, try_strptime(trim(substr(line, 195, 8)), '%Y%m%d')::DATE AS last_observation_date FROM read_csv( getvariable('mpcorb_url'), delim = '|', header = false, columns = {'line': 'VARCHAR'}, quote = '', escape = '', strict_mode = false ) WHERE length(line) >= 202 ), observation_spans AS ( SELECT *, CASE WHEN regexp_full_match(observation_arc, '[0-9]{4}-[0-9]{4}') THEN try_cast(right(observation_arc, 4) AS INTEGER) - try_cast(left(observation_arc, 4) AS INTEGER) ELSE 0 END AS observation_years FROM source_orbits ), eligible_orbits AS ( SELECT * FROM observation_spans WHERE opposition_count >= getvariable('min_oppositions') AND observation_years >= getvariable('min_observation_years') AND eccentricity >= 0.0 AND eccentricity < 1.0 AND semi_major_axis_au > 0.0 AND mean_anomaly_deg IS NOT NULL AND mean_daily_motion_deg IS NOT NULL AND argument_of_perihelion_deg IS NOT NULL AND ascending_node_deg IS NOT NULL AND inclination_deg IS NOT NULL ), epoch_decoded AS ( SELECT *, 100 * (ascii(substr(packed_epoch, 1, 1)) - ascii('A') + 10) + try_cast(substr(packed_epoch, 2, 2) AS INTEGER) AS epoch_year, packed_value(substr(packed_epoch, 4, 1)) AS epoch_month, packed_value(substr(packed_epoch, 5, 1)) AS epoch_day FROM eligible_orbits ), elements AS ( SELECT orbit.*, earth.observation_utc, earth.earth_x_au, earth.earth_y_au, earth.earth_z_au, radians(argument_of_perihelion_deg) AS argument_of_perihelion_rad, radians(ascending_node_deg) AS ascending_node_rad, radians(inclination_deg) AS inclination_rad, radians( mod( mod( mean_anomaly_deg + mean_daily_motion_deg * ( earth.jd_ut + earth.tdb_minus_ut_seconds / 86400.0 - (julian(make_date(epoch_year, epoch_month, epoch_day)) - 0.5) ), 360.0 ) + 360.0, 360.0 ) ) AS mean_anomaly_rad FROM epoch_decoded AS orbit CROSS JOIN earth_ephemeris AS earth ), kepler_solution AS ( SELECT *, list_reduce( range(12)::DOUBLE[], lambda eccentric_anomaly, ignored: eccentric_anomaly - ( eccentric_anomaly - eccentricity * sin(eccentric_anomaly) - mean_anomaly_rad ) / (1.0 - eccentricity * cos(eccentric_anomaly)), CASE WHEN eccentricity < 0.8 THEN mean_anomaly_rad ELSE pi() END ) AS eccentric_anomaly_rad FROM elements ), orbital_plane AS ( SELECT *, semi_major_axis_au * (cos(eccentric_anomaly_rad) - eccentricity) AS x_orbital_au, semi_major_axis_au * sqrt(1.0 - eccentricity * eccentricity) * sin(eccentric_anomaly_rad) AS y_orbital_au FROM kepler_solution ), heliocentric AS ( SELECT *, ( cos(argument_of_perihelion_rad) * cos(ascending_node_rad) - sin(argument_of_perihelion_rad) * sin(ascending_node_rad) * cos(inclination_rad) ) * x_orbital_au + ( -sin(argument_of_perihelion_rad) * cos(ascending_node_rad) - cos(argument_of_perihelion_rad) * sin(ascending_node_rad) * cos(inclination_rad) ) * y_orbital_au AS asteroid_x_au, ( cos(argument_of_perihelion_rad) * sin(ascending_node_rad) + sin(argument_of_perihelion_rad) * cos(ascending_node_rad) * cos(inclination_rad) ) * x_orbital_au + ( -sin(argument_of_perihelion_rad) * sin(ascending_node_rad) + cos(argument_of_perihelion_rad) * cos(ascending_node_rad) * cos(inclination_rad) ) * y_orbital_au AS asteroid_y_au, sin(argument_of_perihelion_rad) * sin(inclination_rad) * x_orbital_au + cos(argument_of_perihelion_rad) * sin(inclination_rad) * y_orbital_au AS asteroid_z_au FROM orbital_plane ), distances AS ( SELECT *, sqrt( pow(asteroid_x_au - earth_x_au, 2) + pow(asteroid_y_au - earth_y_au, 2) + pow(asteroid_z_au - earth_z_au, 2) ) AS distance_au FROM heliocentric ) SELECT observation_utc, packed_designation, readable_designation, absolute_magnitude_h, observation_count, opposition_count, observation_arc, observation_years, last_observation_date, distance_au, distance_au * 149597870.7 AS distance_km, abs( eccentric_anomaly_rad - eccentricity * sin(eccentric_anomaly_rad) - mean_anomaly_rad ) AS kepler_residual FROM distances;
SELECT * FROM earth_ephemeris;
SELECT observation_utc, readable_designation, round(distance_au, 9) AS distance_au, round(distance_km) AS distance_km, observation_arc, opposition_count, kepler_residual FROM asteroid_distances WHERE packed_designation = '~0K8Q';
SELECT readable_designation, round(distance_au, 9) AS distance_au, round(distance_km) AS distance_km, observation_count, opposition_count, observation_arc, last_observation_date, kepler_residual FROM asteroid_distances ORDER BY distance_au LIMIT 10;
直接在 DuckDB Web Shell 中运行：https://shell.duckdb.org/#queries=v0,LOAD-icu~,SET-VARIABLE-min_oppositions-%3D-2~,SET-VARIABLE-min_observation_years-%3D-2~,SET-VARIABLE-mpcorb_url-%3D-'https%3A%2F%2Fwww.minorplanetcenter.net%2Fiau%2FMPCORB%2FMPCORB.DAT.gz'~,SET-VARIABLE-observation_utc-%3D-date_trunc('second'%2C-current_timestamp-AT-TIME-ZONE-'UTC')~,SET-VARIABLE-encoded_observation_utc-%3D-replace(-strftime(-'%25Y%20%25m%20%25dT%25H%3A%25M%3A%25S'%2C-getvariable('observation_utc')-)%2C-'%3A'%2C-'%253A'-)~,SET-VARIABLE-miriade_url-%3D-'https%3A%2F%2Fssp.imcce.fr%2Fwebservices%2Fmiriade%2Fapi%2Fephemcc.php%3F%20name%3Dp%253AEarth%26%20ep%3D'-%7C%7C-getvariable('encoded_observation_utc')-%7C%7C-'%26%20nbd%3D1%26%20tscale%3DUTC%26%20observer%3D%2540sun%26%20theory%3DINPOP%26%20teph%3D4%26%20tcoor%3D2%26%20rplane%3D2%26%20mime%3Djson%26%20output%3D%20%20jd%252C%20%20ttmutc%26%20from%3Dduckdb'~,CREATE-TEMP-TABLE-earth_ephemeris-AS-WITH-response-AS-(-SELECT-d.Date%3A%3ADOUBLE-AS-jd_ut%2C-d.%22TT%20UTC%22%3A%3ADOUBLE-AS-tt_minus_ut_seconds%2C-d.px%3A%3ADOUBLE-AS-earth_x_au%2C-d.py%3A%3ADOUBLE-AS-earth_y_au%2C-d.pz%3A%3ADOUBLE-AS-earth_z_au-FROM-read_json_auto(getvariable('miriade_url'))-AS-payload%2C-UNNEST(payload.data)-AS-item(d)-)%2C-time_scales-AS-(-SELECT-*%2C-radians(-357.53-%2B-0.9856003-*-(jd_ut-%20-2451545.0)-)-AS-earth_mean_anomaly_rad-FROM-response-)-SELECT-getvariable('observation_utc')-AS-observation_utc%2C-jd_ut%2C-tt_minus_ut_seconds-%2B-0.001657-*-sin(earth_mean_anomaly_rad)-%2B-0.000022-*-sin(2.0-*-earth_mean_anomaly_rad)-AS-tdb_minus_ut_seconds%2C-earth_x_au%2C-earth_y_au%2C-earth_z_au-FROM-time_scales~,CREATE-TEMP-MACRO-packed_value(packed_character)-AS-(-CASE-WHEN-packed_character-BETWEEN-'0'-AND-'9'-THEN-ascii(packed_character)-%20-ascii('0')-ELSE-ascii(packed_character)-%20-ascii('A')-%2B-10-END-)~,CREATE-TEMP-TABLE-asteroid_distances-AS-WITH-source_orbits-AS-(-SELECT-trim(substr(line%2C-1%2C-7))-AS-packed_designation%2C-try_cast(trim(substr(line%2C-9%2C-5))-AS-DOUBLE)-AS-absolute_magnitude_h%2C-trim(substr(line%2C-21%2C-5))-AS-packed_epoch%2C-try_cast(trim(substr(line%2C-27%2C-9))-AS-DOUBLE)-AS-mean_anomaly_deg%2C-try_cast(trim(substr(line%2C-38%2C-9))-AS-DOUBLE)-AS-argument_of_perihelion_deg%2C-try_cast(trim(substr(line%2C-49%2C-9))-AS-DOUBLE)-AS-ascending_node_deg%2C-try_cast(trim(substr(line%2C-60%2C-9))-AS-DOUBLE)-AS-inclination_deg%2C-try_cast(trim(substr(line%2C-71%2C-9))-AS-DOUBLE)-AS-eccentricity%2C-try_cast(trim(substr(line%2C-81%2C-11))-AS-DOUBLE)-AS-mean_daily_motion_deg%2C-try_cast(trim(substr(line%2C-93%2C-11))-AS-DOUBLE)-AS-semi_major_axis_au%2C-try_cast(trim(substr(line%2C-118%2C-5))-AS-INTEGER)-AS-observation_count%2C-try_cast(trim(substr(line%2C-124%2C-3))-AS-INTEGER)-AS-opposition_count%2C-trim(substr(line%2C-128%2C-9))-AS-observation_arc%2C-trim(substr(line%2C-167%2C-28))-AS-readable_designation%2C-try_strptime(trim(substr(line%2C-195%2C-8))%2C-'%25Y%25m%25d')%3A%3ADATE-AS-last_observation_date-FROM-read_csv(-getvariable('mpcorb_url')%2C-delim-%3D-'%7C'%2C-header-%3D-false%2C-columns-%3D-%7B'line'%3A-'VARCHAR'%7D%2C-quote-%3D-''%2C-escape-%3D-''%2C-strict_mode-%3D-false-)-WHERE-length(line)-%3E%3D-202-)%2C-observation_spans-AS-(-SELECT-*%2C-CASE-WHEN-regexp_full_match(observation_arc%2C-'%5B0%209%5D%7B4%7D%20%5B0%209%5D%7B4%7D')-THEN-try_cast(right(observation_arc%2C-4)-AS-INTEGER)-%20-try_cast(left(observation_arc%2C-4)-AS-INTEGER)-ELSE-0-END-AS-observation_years-FROM-source_orbits-)%2C-eligible_orbits-AS-(-SELECT-*-FROM-observation_spans-WHERE-opposition_count-%3E%3D-getvariable('min_oppositions')-AND-observation_years-%3E%3D-getvariable('min_observation_years')-AND-eccentricity-%3E%3D-0.0-AND-eccentricity-%3C-1.0-AND-semi_major_axis_au-%3E-0.0-AND-mean_anomaly_deg-IS-NOT-NULL-AND-mean_daily_motion_deg-IS-NOT-NULL-AND-argument_of_perihelion_deg-IS-NOT-NULL-AND-ascending_node_deg-IS-NOT-NULL-AND-inclination_deg-IS-NOT-NULL-)%2C-epoch_decoded-AS-(-SELECT-*%2C-100-*-(ascii(substr(packed_epoch%2C-1%2C-1))-%20-ascii('A')-%2B-10)-%2B-try_cast(substr(packed_epoch%2C-2%2C-2)-AS-INTEGER)-AS-epoch_year%2C-packed_value(substr(packed_epoch%2C-4%2C-1))-AS-epoch_month%2C-packed_value(substr(packed_epoch%2C-5%2C-1))-AS-epoch_day-FROM-eligible_orbits-)%2C-elements-AS-(-SELECT-orbit.*%2C-earth.observation_utc%2C-earth.earth_x_au%2C-earth.earth_y_au%2C-earth.earth_z_au%2C-radians(argument_of_perihelion_deg)-AS-argument_of_perihelion_rad%2C-radians(ascending_node_deg)-AS-ascending_node_rad%2C-radians(inclination_deg)-AS-inclination_rad%2C-radians(-mod(-mod(-mean_anomaly_deg-%2B-mean_daily_motion_deg-*-(-earth.jd_ut-%2B-earth.tdb_minus_ut_seconds-%2F-86400.0-%20-(julian(make_date(epoch_year%2C-epoch_month%2C-epoch_day))-%20-0.5)-)%2C-360.0-)-%2B-360.0%2C-360.0-)-)-AS-mean_anomaly_rad-FROM-epoch_decoded-AS-orbit-CROSS-JOIN-earth_ephemeris-AS-earth-)%2C-kepler_solution-AS-(-SELECT-*%2C-list_reduce(-range(12)%3A%3ADOUBLE%5B%5D%2C-lambda-eccentric_anomaly%2C-ignored%3A-eccentric_anomaly-%20-(-eccentric_anomaly-%20-eccentricity-*-sin(eccentric_anomaly)-%20-mean_anomaly_rad-)-%2F-(1.0-%20-eccentricity-*-cos(eccentric_anomaly))%2C-CASE-WHEN-eccentricity-%3C-0.8-THEN-mean_anomaly_rad-ELSE-pi()-END-)-AS-eccentric_anomaly_rad-FROM-elements-)%2C-orbital_plane-AS-(-SELECT-*%2C-semi_major_axis_au-*-(cos(eccentric_anomaly_rad)-%20-eccentricity)-AS-x_orbital_au%2C-semi_major_axis_au-*-sqrt(1.0-%20-eccentricity-*-eccentricity)-*-sin(eccentric_anomaly_rad)-AS-y_orbital_au-FROM-kepler_solution-)%2C-heliocentric-AS-(-SELECT-*%2C-(-cos(argument_of_perihelion_rad)-*-cos(ascending_node_rad)-%20-sin(argument_of_perihelion_rad)-*-sin(ascending_node_rad)-*-cos(inclination_rad)-)-*-x_orbital_au-%2B-(-%20sin(argument_of_perihelion_rad)-*-cos(ascending_node_rad)-%20-cos(argument_of_perihelion_rad)-*-sin(ascending_node_rad)-*-cos(inclination_rad)-)-*-y_orbital_au-AS-asteroid_x_au%2C-(-cos(argument_of_perihelion_rad)-*-sin(ascending_node_rad)-%2B-sin(argument_of_perihelion_rad)-*-cos(ascending_node_rad)-*-cos(inclination_rad)-)-*-x_orbital_au-%2B-(-%20sin(argument_of_perihelion_rad)-*-sin(ascending_node_rad)-%2B-cos(argument_of_perihelion_rad)-*-cos(ascending_node_rad)-*-cos(inclination_rad)-)-*-y_orbital_au-AS-asteroid_y_au%2C-sin(argument_of_perihelion_rad)-*-sin(inclination_rad)-*-x_orbital_au-%2B-cos(argument_of_perihelion_rad)-*-sin(inclination_rad)-*-y_orbital_au-AS-asteroid_z_au-FROM-orbital_plane-)%2C-distances-AS-(-SELECT-*%2C-sqrt(-pow(asteroid_x_au-%20-earth_x_au%2C-2)-%2B-pow(asteroid_y_au-%20-earth_y_au%2C-2)-%2B-pow(asteroid_z_au-%20-earth_z_au%2C-2)-)-AS-distance_au-FROM-heliocentric-)-SELECT-observation_utc%2C-packed_designation%2C-readable_designation%2C-absolute_magnitude_h%2C-observation_count%2C-opposition_count%2C-observation_arc%2C-observation_years%2C-last_observation_date%2C-distance_au%2C-distance_au-*-149597870.7-AS-distance_km%2C-abs(-eccentric_anomaly_rad-%20-eccentricity-*-sin(eccentric_anomaly_rad)-%20-mean_anomaly_rad-)-AS-kepler_residual-FROM-distances~,SELECT-*-FROM-earth_ephemeris~,SELECT-observation_utc%2C-readable_designation%2C-round(distance_au%2C-9)-AS-distance_au%2C-round(distance_km)-AS-distance_km%2C-observation_arc%2C-opposition_count%2C-kepler_residual-FROM-asteroid_distances-WHERE-packed_designation-%3D-'%3B0K8Q'~,SELECT-readable_designation%2C-round(distance_au%2C-9)-AS-distance_au%2C-round(distance_km)-AS-distance_km%2C-observation_count%2C-opposition_count%2C-observation_arc%2C-last_observation_date%2C-kepler_residual-FROM-asteroid_distances-ORDER-BY-distance_au-LIMIT-10~
为啥识别不了完整链接
