calculates the sun position in a given pos(lat,long) and time. See references inside code

sunPosition.py

# -*- coding: utf-8 -*-
import math 
import datetime

#script for getting the sun position




def sunPosition(
    _year,_month,_day,_hour,
    _min = 0,
    _sec = 0 ,
    _lat = 35.68,
    _long=139.69,
    _timeZone = 9.0
    ):
    '''
    returns the azimuth and elevation of the
    sun in a givent position and date time.

    references:
    *stackoverflow*
    http://stackoverflow.com/questions/8708048/position-of-the-sun-given-time-of-day-latitude-and-longitude

    *天体の位置計算*
    http://www.amazon.co.jp/s/ref=nb_sb_noss_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&url=search-alias%3Daps&field-keywords=%E5%A4%A9%E4%BD%93%E3%81%AE%E4%BD%8D%E7%BD%AE%E8%A8%88%E7%AE%97
    '''
    #get day of the year Œ³’U‚©‚ç‚Ì“ú”
    monthDays = [0,31,28,31,30,31,30,31,31,30,31,30]
    day = _day+sum(monthDays[:_month])
    
    #leap year ‚¤‚邤”N
    if (_year%4==0) and ((_year%400==0) or (_year%100!=0)) and (day >= 60):
        day += 1
     
    #Julian Date - 2400000
    hour = _hour + (float(_min)/60) + (float(_sec)/3600) - _timeZone
    delta = _year - 1949
    leap = math.trunc(delta/4)
    jd = 32916.5 + delta * 365 + leap + day + hour / 24
    
    time = jd - 51545.0
    
    #Ecliptic coordinates ‰©“¹À•W
    #Mean longitude •½‹Ï‰©Œo
    mnlong = 280.460 + .9856474 * time
    mnlong %= 360.0
    
    #Mean anomaly •½‹Ï‹ß“_Šp
    mnanom = 357.528 + .9856003 * time
    mnanom %= 360.0
    mnanom = math.radians(mnanom)
    
    #Ecliptic longitude and obliquity of ecliptic
    eclong = mnlong+1.915*math.sin(mnanom) + 0.020 * math.sin(2*mnanom)
    eclong %= 360.0
    #‰©Œo: ‰©“¹‚ðŠî€‚Æ‚µ‚½Œo“xAt•ª“_‚ð0‚Æ‚µ‚Ä360‚܂ŁAH‚Í180“x
    eclong = math.radians(eclong)
    
    #angle of equatorial plan and ecpliptic plane ’n‹…‚ÌŒX‚«
    oblqec = 23.439 - 0.0000004 * time 
    oblqec = math.radians(oblqec)
    
    #Celestial coordinates “V‘̍À•W
    #Right ascention and declination ÔŒo ‚Æ ÔˆÜ
    num = math.cos(oblqec) * math.sin(eclong)
    den = math.cos(eclong)
    ra = math.atan(num / den) # <- ÔŒo
    if den < 0 : ra += math.pi
    if (den >= 0) and (num < 0): ra += math.pi*2
    dec = math.asin(math.sin(oblqec)* math.sin(eclong)) # <- ÔˆÜ
    
    #Local coordinates
    #greenwich mean sidereal time
    gmst = (6.697375 + .0657098242 * time + hour) % 24
    #local mean sidereal time
    lmst = (gmst + _long/15.0) % 24
    lmst = math.radians(lmst*15.0)
    
    #hour angle ŽžŠp
    ha = lmst - ra
    if ha < -(math.pi):  ha += math.pi*2
    elif ha > math.pi: ha -= math.pi*2
    
    #latitude to radians
    lat = math.radians(_lat)

    #azimuth and elevation
    el = math.asin(
    math.sin(dec)*math.sin(lat)+
    math.cos(dec)*math.cos(lat)*math.cos(ha)
    ) #‚“xi‹ÂŠpj
    
    az = math.asin(
    -math.cos(dec)*math.sin(ha)/math.cos(el)
    ) #•ûˆÊŠp
    
    cosAzPos = 0 <= math.sin(dec) - math.sin(el) * math.sin(lat)
    sinAzNeg = math.sin(az) < 0
    
    if cosAzPos and sinAzNeg : az += math.pi*2
    if not cosAzPos : az = math.pi - az
    
    elevation = math.degrees(el)
    azimuth = math.degrees(az)
    latitude = math.degrees(lat)
    
    return elevation,azimuth
    


if __name__ == "__main__":

    #testing
    year = 2013
    month = 6
    day = 14
    hour = 1
    min = 0
    sec = 0
    lat = 33.59
    long = 130.314

    sunP = sunPosition(year,month,day,hour,min,sec,lat,long)
    print sunP