BEGIN{ d = julianday(year, month, day) - julianday(2000, 1, 1.5) t = 0 h = 280.466449 + 0.985647360 * d s = 218.316656 + 13.17639647754 * d p = 83.353243 + 0.11140352392 * d n = 125.044556 - 0.05295376276 * d t_hour = 15 s_hour = 0.54901652 h_hour = 0.04106864 p_hour = 0.00464181 s_m2 = 2 * t_hour - 2 * s_hour + 2 * h_hour + 0 * p_hour s_k1 = 1 * t_hour + 0 * s_hour + 1 * h_hour + 0 * p_hour s_s2 = 2 * t_hour + 0 * s_hour + 0 * h_hour + 0 * p_hour s_o1 = 1 * t_hour - 2 * s_hour + 1 * h_hour + 0 * p_hour lon = -140.05 time_zone = -9 h_m2 = 50.9 k_m2 = 152.5 h_k1 = 25.7 k_k1 = 178.9 h_s2 = 24.8 k_s2 = 181.3 h_o1 = 20.0 k_o1 = 160.0 z0 = 120 v_m2_0 = v0_m2(t, s, h, p) + u_m2(n) - 2 * lon + s_m2 * time_zone - k_m2 v_k1_0 = v0_k1(t, s, h, p) + u_k1(n) - 1 * lon + s_k1 * time_zone - k_k1 v_s2_0 = v0_s2(t, s, h, p) + u_s2(n) - 2 * lon + s_s2 * time_zone - k_s2 v_o1_0 = v0_o1(t, s, h, p) + u_o1(n) - 1 * lon + s_o1 * time_zone - k_o1 for(hour = 0; hour < 24; hour++){ eta_m2 = f_m2(n) * h_m2 * cos(rad(s_m2 * hour + v_m2_0)) eta_k1 = f_k1(n) * h_k1 * cos(rad(s_k1 * hour + v_k1_0)) eta_s2 = f_s2(n) * h_s2 * cos(rad(s_s2 * hour + v_s2_0)) eta_o1 = f_o1(n) * h_o1 * cos(rad(s_o1 * hour + v_o1_0)) printf "%d-%d-%dT%02d:00:00 %f\n", year, month, day, hour, eta_m2 + eta_k1 + eta_s2 + eta_o1 + z0 } } function v0_m2(t, s, h, p){ return 2 * t - 2 * s + 2 * h + 0 * p + 0 } function v0_k1(t, s, h, p){ return 1 * t + 0 * s + 1 * h + 0 * p + 90 } function v0_s2(t, s, h, p){ return 2 * t + 0 * s + 0 * h + 0 * p + 0 } function v0_o1(t, s, h, p){ return 1 * t - 2 * s + 1 * h + 0 * p + 270 } function f_m2(n){ tmp = 1.0004 tmp += -0.0373 * cos (rad(n)) tmp += 0.0002 * cos(rad(n) * 2) tmp += 0.0000 * cos(rad(n) * 3) return tmp } function u_m2(n){ tmp = -2.14 * sin(rad(n)) tmp += 0.00 * sin(rad(n) * 2) tmp += 0.00 * sin(rad(n) * 3) return tmp } function f_k1(n){ tmp = 1.0060 tmp += 0.1150 * cos(rad(n)) tmp += -0.0088 * cos(rad(n) * 2) tmp += 0.0006 * cos(rad(n) * 3) return tmp } function u_k1(n){ tmp = -8.86 * sin(rad(n)) tmp += 0.68 * sin(rad(n) * 2) tmp += -0.07 * sin(rad(n) * 3) return tmp } function f_s2(n){ return 1 } function u_s2(n){ return 0 } function f_o1(n){ tmp = 1.0089 tmp += 0.1871 * cos(rad(n)) tmp += -0.0147 * cos(rad(n) * 2) tmp += 0.0014 * cos(rad(n) * 3) return tmp } function u_o1(n){ tmp = 10.80 * sin(rad(n)) tmp += -1.34 * sin(rad(n) * 2) tmp += 0.19 * sin(rad(n) * 3) return tmp } function rad(deg){ return deg * atan2(0, -1) / 180 } function julianday(y, m, d){ if(m <= 2){ y -= 1 m += 12 } a = int(y / 100) b = 2 - a + int(a / 4) c = int(365.25 * y) return 1720994 + int(30.6001 * (m + 1)) + b + c + d + 0.5 }