太阳直射点纬度与日期公式，python程序

太阳赤纬=0.3723＋23.2567sinθ＋0.1149sin2θ－0.1712sin3θ－0.758cosθ＋0.3656cos
2θ＋0.0201cos3θ
式中θ称日角，即 θ=2πt／365.2422
这里t又由两部分组成，即 t=N－N0
式中N为积日，所谓积日，就是日期在年内的顺序号，例如，1月1日其积日为1，平年12月
31日的积日为365，闰年则为366，等等。
N0=79.6764＋0.2422×（年份－1985）－INT〔（年份－1985）／4〕

#!/usr/bin/python
import math,os,sys

def jiri(year,month,day):
    if year%4==0 and year%100!=0 or year%400==0:leap=1
    else:leap=0
    mm=(0,31,leap+59,leap+90,leap+120,leap+151,leap+181,leap+212,leap+243,leap+273,leap+304,leap+334)
    return mm[month-1]+day

if __name__=="__main__":
    year=int(sys.argv[1])
    month=int(sys.argv[2])
    day=int(sys.argv[3])
    t=jiri(year,month,day)-79.6764-0.2422*(year-1985)+int((year-1985)/4)
    sita=2*3.14159265*t/365.2422
    delta=0.3723+23.2567*math.sin(sita)+0.1149*math.sin(2*sita)-0.1712*math.sin(3*sita)-0.758*math.cos(sita)+0.3656*math.cos(
2*sita)+0.0201*math.cos(3*sita)
    if delta<0:s='-'
    else:s='+'
    delta=abs(delta)
    d=int(delta)
    f=int((delta-d)*60)
    m=delta*3600-f*60-d*3600
    print s,d,"degree",f,"minute",m,"second"