1       REM
2       REM     Handling program for SUNRS
3       REM
paluu etusivulle
5       DEF FNRAD (W) = .01745329252# * W
10      DEF FNDEG (W) = 57.29577951# * W

15      DIM FLAG(9)
19      CLS
20      PRINT : FLAG(0) = 0
25      PRINT "Sunrise and Sunset V1.0 M.Vanhanen"
30      PRINT "======================================================"

35      IF FLAG(9) = 1 THEN : GOTO 65
37      PRINT "Maantieteellinen pituus, etäisyys nollameridiaanista itään tai länteen"
38      PRINT "Anna negatiivinen arvo (-), jos on länteen meridiaanista"
40      INPUT " (HELSINKI 24ø20'0'') Geog. longitude (W neg; D,M,S)"; XD, XM, XS
41      PRINT
45      SW1 = -1: GOSUB 1000: LNG = X
47      PRINT "Leveysaste, päiväntasaajasta pohjoiseen tai etelään"
48      PRINT "Anna negatiivinen øarvo (-), jos on etelään päiväntasaajasta"
50      INPUT " (HELSINKI 60ø10'0'') Geographical latitude (E neg; D,M,S)"; XD, XM, XS
55      GOSUB 1000: PHI = FNRAD(X)
60      PRINT

65      IF FLAG(1) = 1 THEN : GOTO 80
70      INPUT "Date (D,M,Y) [dd,mm,yyyy]"; dy, mn, yr
71      LOCATE 12, 1: PRINT "Date "; dy; ", Month "; mn; ", Year "; yr; ".           "
75      PRINT

80      GOSUB 3500
85      IF ERR1 <> 0 THEN GOTO 130
90      PRINT "     GMT of sunrise: "; HSR; "h"; MSR; "m"; SSR; "s"
95      IF ERR2 = 1 THEN PRINT "** ambiguous conversion **"
100     PRINT "     ... and sunset: "; HSS; "h"; MSS; "m"; SSS; "s"
102     PRINT : PRINT "Keskiaurinkoaikaa!, Suomessa oltaessa lisää +2h (kesä aikana +3h)."
105     IF ERR3 = 1 THEN PRINT "** ambiguous conversion **"
110     AZR0 = FNDEG(AZR0): AZS0 = FNDEG(AZS0)
115     PRINT : PRINT "  Azimuths (deg) are"
120     PRINT "         ....rising: "; AZR0; "ø"
125     PRINT "    ... and setting: "; AZS0; "ø"
130     PRINT

135     INPUT "Again (Y/N)"; an$
140     IF an$ = "N" OR an$ = "n" THEN STOP
145     INPUT "New place (Y/N)"; an$
150     IF an$ = "Y" OR an$ = "y" THEN FLAG(9) = 0
155     INPUT "New date (Y/N)"; an$
160     IF an$ = "Y" OR an$ = "y" THEN FLAG(1) = 0: FLAG(5) = 0: FLAG(6) = 0: FLAG(3) = 0
165     CLS : GOTO 30

997     REM
998     REM     Subroutine MINSEC
999     REM

1000    IF SW1 = -1 THEN GOTO 1035
1005    SN = SGN(X): XP = ABS(X): XD = INT(XP)
1010    A = (XP - XD) * 60: XM = INT(A)
1015    XS = INT((A - XM) * 600 + .5) / 10
1020    S$ = "+"
1025    IF SN = -1 THEN S$ = "-"
1030    RETURN

1035    SN = 1
1040    IF XD < 0 OR XM < 0 OR XS < 0 THEN SN = -1
1045    XD1 = ABS(XD): XM1 = ABS(XM): XS1 = ABS(XS)
1050    X = ((((XS1 / 60) + XM1) / 60) + XD1) * SN
1055    RETURN

1097    REM
1098    REM     Subroutine JULDAY
1099    REM

1100    IF FLAG(1) = 1 THEN RETURN
1105    DEF FNITG (W) = INT(W) + FLAG(0) * SGN(W) * INT((SGN(W) - 1) / 2) * SGN(W - INT(W))
1110    MN1 = mn: YR1 = yr: FLAG(1) = 1: B = 0
1115    IF YR1 < 0 THEN YR1 = YR1 + 1
1120    IF mn < 3 THEN MN1 = mn + 12: YR1 = YR1 - 1
1125    IF yr < 1582 THEN GOTO 1145
1130    IF yr = 1582 AND mn < 10 THEN GOTO 1145
1135    IF yr = 1582 AND mn = 10 AND dy < 15 THEN GOTO 1145
1140    A = INT(YR1 / 100): B = 2 - A + INT(A / 4)
1145    IF YR1 < 0 THEN GOTO 1155
1150    C = INT(365.25 * YR1) - 694025: GOTO 1160
1155    C = FNITG((365.25 * YR1) - .75) - 694025
1160    D = INT(30.6001 * (MN1 + 1))
1165    DJD = B + C + D + dy - .5
1170    RETURN

1297    REM
1298    REM     Subroutine GTIME
1299    REM

1300    IF FLAG(3) = 1 THEN GOTO 1340
1305    XMN = mn: XDY = dy: FLAG(3) = 1
1310    FLAG(1) = 0: GOSUB 1100: XDJD = DJD
1315    mn = 1: dy = 0: FLAG(1) = 0
1320    GOSUB 1100: T = DJD / 36525
1325    R = 6.6460656# + (.051262 + (T * 2.581E-05)) * T
1330    R1 = 2400 * (T - ((yr - 1900) / 100)): B = 24 - R - R1
1335    T0 = ((XDJD - DJD) * .0657098) - B
1340    IF SW2 = 1 THEN GOTO 1365

1345    IF T0 < 0 THEN T0 = T0 + 24
1350    GTM = TIM - T0
1355    IF GTM < 0 THEN GTM = GTM + 24
1360    GTM = GTM * .9972695677#: GOTO 1380

1365    GTM = (TIM * 1.002737908#) + T0
1370    IF GTM > 24 THEN GTM = GTM - 24
1375    IF GTM < 0 THEN GTM = GTM + 24

1380    X = GTM: SW1 = 1: GOSUB 1000
1385    HS = XD: MS = XM: SS = XS
1390    mn = XMN: dy = XDY: DJD = XDJD
1395    RETURN

1697    REM
1698    REM     Subroutine OBLIQ
1699    REM

1700    IF FLAG(5) = 1 THEN RETURN
1705    FLAG(5) = 1
1710    GOSUB 1100: T = DJD / 36525
1715    C = (((-.00181 * T) + .0059) * T + 46.845) * T
1720    EPS = 23.45229444# - (C / 3600)
1725    EPSR = EPS * .01745329252#
1730    RETURN

1797    REM
1798    REM     Subroutine EQECL
1799    REM

1800    IF FLAG(7) = 1 THEN GOTO 1830
1805    DEF FNASN (W) = ATN(W / (SQR(1 - W * W) + 1E-20))
1810    PI = 3.1415926535#: TPI = 2 * PI: FLAG(7) = 1
1815    IF FLAG(6) = 0 THEN DEPS = 0
1820    GOSUB 1700: EPS1 = FNRAD(EPS + DEPS)
1825    SEPS = SIN(EPS1): CEPS = COS(EPS1)

1830    CY = COS(Y): SY = SIN(Y)
1835    IF ABS(CY) < 1E-20 THEN CY = 1E-20
1840    TY = SY / CY: CX = COS(X): SX = SIN(X)
1845    SQ = (SY * CEPS) - (CY * SEPS * SX * SW3)
1850    Q = FNASN(SQ): A = (SX * CEPS) + (TY * SEPS * SW3)
1855    P = ATN(A / CX)
1860    IF CX < 0 THEN P = P + PI
1865    IF P > TPI THEN P = P - TPI
1870    IF P < 0 THEN P = P + TPI
1875    RETURN

2097    REM
2098    REM     Subroutine RISET
2099    REM

2100    IF FLAG(9) = 1 THEN GOTO 2125
2105    REM DEF FNASN (W) = ATN(W / (SQR(1 - W * W) + 1E-20))
2110    DEF FNACS (W) = 1.570796327# - FNASN(W)
2115    CPHI = COS(PHI): SPHI = SIN(PHI)
2120    TPI = 6.283185308#: FLAG(9) = 1

2125    CY = COS(Y): A = CY * CPHI
2130    IF ABS(A) < 1E-10 THEN GOTO 2260

2135    CPSI = SPHI / CY
2140    IF CPSI > 1 THEN CPSI = 1
2145    IF CPSI < -1 THEN CPSI = -1
2150    PSI = FNACS(CPSI): SPSI = SIN(PSI)
2155    DH = DIS * SPSI: Y1 = Y + (DIS * CPSI)
2160    SY = SIN(Y1): TY = TAN(Y1)
2165    CH = (-1 * SPHI * TY) / CPHI

2170    IF CH < -1 THEN GOTO 2240
2175    IF CH > 1 THEN GOTO 2250

2180    ERR1 = 0: H = FNACS(CH)
2185    H = (H + DH) * 3.819718634#
2190    LSTR = 24 + X - H: LSTS = X + H
2195    CAZR = SY / CPHI: AZR = FNACS(CAZR)
2200    AZS = TPI - AZR: B = LSTR: C = 24: GOSUB 2225
2205    LSTR = B: B = LSTS: GOSUB 2225: LSTS = B
2210    B = AZR: C = TPI: GOSUB 2225: AZR = B
2215    B = AZS: GOSUB 2225: AZS = B
2220    RETURN

2225    IF B < 0 THEN B = B + C
2230    IF B > C THEN B = B - C
2235    RETURN

2240    ERR1 = -1: PRINT "** circumpolar **"
2245    RETURN
2250    ERR1 = 1: PRINT "** never rises **"
2255    RETURN
2260    ERR1 = 2: PRINT "** impossible calculation **"
2265    RETURN

2697    REM
2698    REM     Subroutine NUTAT
2699    REM

2700    IF FLAG(6) = 1 THEN RETURN
2705    REM DEF FNRAD (W) = .01745329252# * W
2710    FLAG(6) = 1: FLAG(7) = 0
2715    GOSUB 1100: T = DJD / 36525: T2 = T * T
2720    A = 100.0021358# * T: B = 360 * (A - INT(A))
2725    LS = 279.6967 + .000303 * T2 + B
2730    A = 1336.855231# * T: B = 360 * (A - INT(A))
2735    LD = 270.4342 - .001133 * T2 + B
2740    A = 99.99736056# * T: B = 360 * (A - INT(A))
2745    MS = 358.4758 - .00015 * T2 + B
2750    A = 1325.552359# * T: B = 360 * (A - INT(A))
2755    MD = 296.1046 + .009192 * T2 + B
2760    A = 5.372616667# * T: B = 360 * (A - INT(A))
2765    NM = 259.1833 + .002078 * T2 - B
2770    TLS = 2 * FNRAD(LS): NM = FNRAD(NM)
2775    TNM = 2 * FNRAD(NM): MS = FNRAD(MS)
2780    TLD = 2 * FNRAD(LD): MD = FNRAD(MD)

2781    DPSI = (-17.2327 - .01737 * T) * SIN(NM)
2782    DPSI = DPSI + (-1.2729 - .00013 * T) * SIN(TLS)
2783    DPSI = DPSI + .2088 * SIN(TNM) - .2037 * SIN(TLD)
2784    DPSI = DPSI + (.1261 - .00031 * T) * SIN(MS)
2785    DPSI = DPSI + .0675 * SIN(MD)
2786    DPSI = DPSI - (.0497 - .00012 * T) * SIN(TLS + MS)
2787    DPSI = DPSI - .0342 * SIN(TLD - NM) - .0261 * SIN(TLD + MD)
2788    DPSI = DPSI + .0214 * SIN(TLS - MS)
2789    DPSI = DPSI - .0149 * SIN(TLS - TLD + MD)
2790    DPSI = DPSI + .0124 * SIN(TLS - NM) + .0114 * SIN(TLD - MD)

2791    DEPS = (9.21 + .00091 * T) * COS(NM)
2792    DEPS = DEPS + (.5522 - .00029 * T) * COS(TLS)
2793    DEPS = DEPS - .0904 * COS(TNM) + .0884 * COS(TLD)
2794    DEPS = DEPS + .0216 * COS(TLS + MS) + .0183 * COS(TLD - NM)
2795    DEPS = DEPS + .0113 * COS(TLD + MD) - .0093 * COS(TLS - MS)
2796    DEPS = DEPS - .0066 * COS(TLS - NM)

2797    DPSI = DPSI / 3600: DEPS = DEPS / 3600
2800    RETURN

3097    REM
3098    REM     Subroutine ANOMALY
3099    REM

3100    TPI = 6.283185308#
3105    M = MA - TPI * INT(MA / TPI): EA = M
3110    DLA = EA - (S * SIN(EA)) - M
3115    IF ABS(DLA) < .000001 THEN GOTO 3130
3120    DLA = DLA / (1 - (S * COS(EA)))
3125    EA = EA - DLA: GOTO 3110

3130    TNU2 = SQR((1 + S) / (1 - S)) * TAN(EA / 2)
3135    NU = 2 * ATN(TNU2)
3140    RETURN

3297    REM
3298    REM     Subroutine SUN
3299    REM

3300    GOSUB 1100: T = DJD / 36525: T2 = T * T
3305    REM DEF FNRAD (W) = .01745329252# * W
3310    A = 100.0021359# * T: B = 360 * (A - INT(A))
3315    LS = 279.69668# + .0003025 * T2 + B
3320    A = 99.99736042000001# * T: B = 360 * (A - INT(A))
3325    MS = 358.47583# - (.00015 + .0000033 * T) * T2 + B
3330    S = 1.675104E-02 - .0000418 * T - 1.26E-07 * T2
3335    MA = FNRAD(MS): GOSUB 3100
3340    A = 62.55209472# * T: B = 360 * (A - INT(A))
3345    A1 = FNRAD(153.23 + B)
3350    A = 125.1041894# * T: B = 360 * (A - INT(A))
3355    B1 = FNRAD(216.57 + B)
3360    A = 91.56766028# * T: B = 360 * (A - INT(A))
3365    C1 = FNRAD(312.69 + B)
3370    A = 1236.853095# * T: B = 360 * (A - INT(A))
3375    D1 = FNRAD(350.74 - .00144 * T2 + B)
3380    E1 = FNRAD(231.19 + 20.2 * T)
3385    A = 183.1353208# * T: B = 360 * (A - INT(A))
3390    H1 = FNRAD(353.4 + B)
3395    DL = .00134 * COS(A1) + .00154 * COS(B1)
3396    DL = DL + .002 * COS(C1) + .00179 * SIN(D1)
3397    DL = DL + .00178 * SIN(E1)
3400    DR = 5.43E-06 * SIN(A1) + 1.575E-05 * SIN(B1)
3401    DR = DR + 1.627E-05 * SIN(C1) + 3.076E-05 * COS(D1)
3402    DR = DR + 9.27E-06 * SIN(H1)
3405    LSN = NU + FNRAD(LS - MS + DL): TPI = 6.283185308#
3410    RSN = 1.0000002# * (1 - S * COS(EA)) + DR
3415    IF LSN < 0 THEN LSN = LSN + TPI: GOTO 3415
3420    IF LSN > TPI THEN LSN = LSN - TPI: GOTO 3420
3425    RETURN

3497    REM
3498    REM     Subroutine SUNRS
3499    REM

3500    DY0 = dy: dy = DY0 + .5
3505    DIS = .014593: GOSUB 3645
3510    IF ERR1 <> 0 THEN RETURN
3515    AL = LNG / 15: TR = LSTR - AL: TS = LSTS - AL
3520    GOSUB 3620: dy = DY0: TIM = TR: SW2 = -1
3525    GOSUB 1300: ERR2 = 0: ERR3 = 0
3530    IF GTM < .06552 THEN ERR2 = 1
3535    GTMR = GTM: TIM = TS: GOSUB 1300
3540    IF GTM < .06552 THEN ERR3 = 1
3545    GTMS = GTM: DJD0 = DJD
3550    DJD = DJD0 + GTMR / 24: GOSUB 3645
3555    IF ERR1 <> 0 THEN RETURN
3560    AZR0 = AZR: TR = LSTR - AL
3565    DJD = DJD0 + GTMS / 24: GOSUB 3645
3570    IF ERR1 <> 0 THEN RETURN
3580    AZS0 = AZS: TS = LSTS - AL: GOSUB 3620
3585    TIM = TR: GOSUB 1300
3590    IF GTM < .06552 THEN ERR2 = 1
3595    HSR = HS: MSR = MS: SSR = SS
3600    TIM = TS: GOSUB 1300
3605    IF GTM < .06552 THEN ERR3 = 1
3610    HSS = HS: MSS = MS: SSS = SS
3615    RETURN

3620    IF TR < 0 THEN TR = TR + 24
3625    IF TS < 0 THEN TS = TS + 24
3630    IF TR > 24 THEN TR = TR - 24
3635    IF TS > 24 THEN TS = TS - 24
3640    RETURN

3645    GOSUB 3300: GOSUB 2700: HT = 0
3650    X = LSN + .01745329252# * (DPSI - .00569)
3655    Y = 0: SW3 = -1: GOSUB 1800
3660    X = 3.819718634# * P: Y = Q: GOSUB 2100
3665    RETURN
paluu etusivulle