!日の出、日の入りの計算（Sunrise equation）
!西暦2010年11月16日
!東経139.7447、北緯35.6544
!日の出06:17:06、方位角112.481
!日の入16:34:01、方位角247.3833
OPTION ANGLE DEGREES
LET year=2010 !グレゴリオ暦
LET month=11 !月
LET day=16 !日
LET lw=-139.7447 !経度（単位は度） ※東経は負
LET phi=35.6544 !緯度（単位は度） ※北緯は正
!Converting Gregorian calendar date to Julian Day Number
LET Jd=JDN(year,month,day)
!CALL TransDateToJD(year,month,day,HH,MM,SS, jd)
PRINT Jd !debug
!Calculate current Julian Cycle（2000/1/1=2451545）
PRINT JDN(2000,1,1) !debug
LET nn=Jd-JDN(2000,1,1)-0.0009-lw/360
LET n=ROUND(nn,0)
!Approximate Solar Noon
LET J=JDN(2000,1,1)+0.0009+lw/360+n
!Solar Mean Anomaly
LET M=MOD( 357.5291+0.98560028*(J-JDN(2000,1,1)) , 360 )
!Equation of Center
LET C=1.9148*SIN(M)+0.0200*SIN(2*M)+0.0003*SIN(3*M)
!Ecliptic Longitude
LET lambda=MOD( M+102.9372+C+180 , 360 )
!Solar Transit
LET Jt=J+0.0053*SIN(M)-0.0069*SIN(2*lambda)
!Declination of the Sun（太陽赤緯）
LET delta=ASIN( SIN(lambda)*SIN(23.45) )
!Hour Angle（時角）
LET w0=ACOS( (SIN(-0.83)-SIN(phi)*SIN(delta))/(COS(phi)*COS(delta)) )
!Calculate Sunrise and Sunset
LET Jset=JDN(2000,1,1)+0.0009+((w0+lw)/360+n+0.0053*SIN(M))-0.0069*SIN(2*lambda)
LET Jrise=Jt-(Jset-Jt)
PRINT Jset; Jrise !debug
CALL TransJDToDate(Jrise, Y,M,D,HH,MM,SS)
PRINT Y;"年";M;"月";D;"日"; HH;"時";MM;"分";SS;"秒"
CALL TransJDToDate(Jset, Y,M,D,HH,MM,SS)
PRINT Y;"年";M;"月";D;"日"; HH;"時";MM;"分";SS;"秒"
END
EXTERNAL FUNCTION JDN(year,month,day) !Converting Gregorian calendar date to Julian Day Number
LET a=INT((14-month)/12)
LET y=year+4800-a
LET m=month+12*a-3
LET JDN=day+INT((153*m+2)/5)+365*y+INT(y/4)-INT(y/100)+INT(y/400)-32045
END FUNCTION
EXTERNAL SUB TransDateToJD(Y,M,D,HH,MM,SS, jd) !グレゴリオ暦（日時）からユリウス日を得る
IF M<3 THEN
LET Y=Y-1
LET M=M+12
END IF
LET t =INT(Y/100)
LET jd=INT(Y*365.25)-t+INT(t/4)
LET jd=jd+INT(30.6001*(M+1))+D+HH/24+MM/1440+SS/86400+1720996.5
LET jd=jd-9/24 !JST、UTならNOP
END SUB
EXTERNAL SUB TransJDToDate(jd, Y,M,D,HH,MM,SS) !ユリウス日からグレゴリオ暦（日時）を得る
LET jd=jd+0.875 !JST、UTなら0.5
LET z =INT(jd)
LET f =jd-z
LET aa=INT((z-1867216.25)/36524.25)
LET a =INT(z+1+aa-INT(aa/4))
LET b =a+1524
LET c =INT((b-122.1)/365.25)
LET k =INT(365.25*c)
LET e =INT((b-k)/30.6001)
LET D=INT(b-k-INT(30.6001*e))
IF e<13.5 THEN LET M=e-1 ELSE LET M=e-13
IF M>2.5 THEN LET Y=c-4716 ELSE LET Y=c-4715
LET HH=INT(f*24)
LET MM=INT((f*24-HH)*60)
LET SS=((f*24-HH)*60-MM)*60
END SUB