# coding: utf-8
# ***** BEGIN GPL LICENSE BLOCK *****
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#
# ***** END GPL LICENCE BLOCK *****
# Planetarium Unit 2008.11.25(Tue) 01:47
import sys,re,math,time
from math import pi,cos,sin,atan,acos,asin,tan,floor,sqrt,atan2

_n_=0 ; _s_=1 ; _e_=2 ; _w_=3
_NSEW_={'n':_n_ , 'N':_n_,
        's':_s_ , 'S':_s_,
	'e':_e_ , 'E':_e_,
	'w':_w_ , 'W':_w_ }
_rexp_=re.compile("[^\s]+")

_iDir_=[1,-1]

def Sind(deg) : return sin(deg*(pi/180.0))
def Cosd(deg) : return cos(deg*(pi/180.0))
def aSind(y)  : return asin(y)/(pi/180.0)
def aCosd(y)  : return acos(y)/(pi/180.0)

# + a Julian Day from a time tuple
def JulianDay(year=0, month=1, day=1, hour=0, minute=0, second=0):
  # Algorithm:
  # Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,
  # Virginia. p. 61
  # Convert time to fractions of a day
  day = day + hour/24.0 + minute/1440.0 + second/86400.0 
  
  # Start Meeus algorithm (variables are in his notation)
  if (month < 3):
    month = month + 12
    year = year -1
  
  A = int(year/100.0)
  jd = int(365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + day - 1524.5
  
  # Adjust the jd for the switch from the Julian to Gregorian Calendar
  # here assumed to have occurred the day after 1582 October 4
  if jd > 2299169.5:
    # 1582 October 15 (Gregorian Calendar)
    B = 2 - A + int(A/4)
  elif jd < 2299160.5:
    # 1582 October 4 (Julian Calendar)
    B = 0
  else:
    date = "%04d-%02d-%02d %02d:%02d:%02d" % (year, month, day, hour, minute, second)
  # adjust for Julian calendar if necessary
  jd = jd + B
  return jd

# + Sidereal Time
def GST(y,m,d,hr,mn,sc):
  j=JulianDay(y,m,d,hr,mn,sc)
  T=(j-2451545)/36525.0
  G0=24110.548 + (8640184.812866*T)+(0.93104*T*T)-(0.0000062*T*T*T)
  a=(hr*3600+mn*60+sc)*1.00273790934
  return (((a+G0)/3600.0)+9600)%24.0

# Determine the Local Sidereal Date and Time from the longitude
def LST(y,m,d,hr,mn,sc,L):
  return GST(y,m,d,hr,mn,sc)-(L/15.0)

_sxg_=0 ; _hra_=1 ; _pct_=2 ; _rad_=3 ;
_rdAngle_=[1,15,3.6,(180/pi)]

# 2007.08.09(Thu) 14:41 Decimal process
def rdAngle(v,flag=_sxg_):
  v = _rexp_.findall(v)
  nsew=-1
  if len(v)>3:
    # n:0 s:1 e:2 w=3
    if _NSEW_.has_key(v[3]):
      nsew=_NSEW_[v[3]]
  # str to int
  dd=int(v[0])
  mm=int(v[1])
  ss=int(float(v[2]))
  pm=_iDir_[int(dd)<0]
  
  # + decimal 60 to 10
  #   float problem
  ang=abs(dd)+(abs(mm)/60.0)+(abs(ss)/3600.0)
  ang=pm*ang
  
  # + W N convert
  #   - west longitude
  if (nsew==_e_) : ang=360-ang
  #   - southern latitude
  if (nsew==_s_) : ang=0-ang
  
  ang*=_rdAngle_[flag]
  return ang

# + Equatorial coordinate -> coordinate the horizons of observation points
#   Greenwich mean sidereal time(2000.0)[deg]
#   JD:Julian Day
#   alpha:Right ascension[rad]  delta:Declination[rad]
#   Lo(longitude):East [deg]
#   la(latitude) :North[deg]
#   az:Azimuth[deg]  el:Elevation[deg]
def horizontal(TM,JD,alpha,delta,Lo,la,hr,mn,sc) :
  # 1 1900.1.0 noon  days elapsed
  K    = JD - 2415020.0
  # 2 100 years / unit
  Tu   = K/36525.0
  # 3 GMST=Greenwich mean sidereal time
  gmst = 6+(38.0/60.0)+(45.836/3600.0) + 8640184.542*Tu/3600 + 0.0929*pow(Tu,2)/3600
  while(gmst>24): gmst-=24
  # - deg convert
  gmst *= 15
  # + Apparent sidereal time[deg]
  TIME_DIFFERENCE=5.5
  hr = 15*1.00273791*( hr+mn/60.0+sc/3600.0 - TIME_DIFFERENCE )
  st = gmst + Lo + hr
  #   - range 0..360
  while (st>360):st-=360
  while (st<=0) :st+=360
  
  # 5 alpha?蝠景ght ascension delta?蜥ｼeclination
  # + to rad
  st *= pi/180.0
  la *= pi/180.0
  # + Horizontal coordinates
  x =  sin(la)*sin(delta) + cos(la)*cos(delta)*cos(st-alpha)
  y = -cos(la)*sin(delta) + sin(la)*cos(delta)*cos(st-alpha)
  z =                               cos(delta)*sin(st-alpha)
  # Elevation
  cos_el = sqrt( 1 - pow(x,2) )
  el     = acos( cos_el)			# rad to deg
  el     *= 180/pi				#
  # Azimuth
  az     = atan2(z,y)
  az     *= 180/pi
  if( az<0 ) : az += 360.0
  
  return x,y,z,az,el

# + Star Color and Size
_StarSpec_={
'O': [ float(155)/255, float(176)/255, float(255)/255 ],
'B': [ float(170)/255, float(191)/255, float(255)/255 ],
'A': [ float(202)/255, float(215)/255, float(255)/255 ],
'F': [ float(248)/255, float(247)/255, float(255)/255 ],
'G': [ float(255)/255, float(244)/255, float(234)/255 ],
'K': [ float(255)/255, float(210)/255, float(161)/255 ],
'M': [ float(255)/255, float(204)/255, float(111)/255 ],
'N': [ float(232)/255,   0.0,   0.0 ],
'C': [ float(210)/255,   0.0,   0.0 ],
'S': [ float(195)/255,   0.0,   0.0 ],
'L': [ float(128)/255,   0.0,   0.0 ],
'T': [ float( 32)/255,   0.0,   0.0 ],
'-': [ 0.0,0.0,0.0]
}

_StarSize_=[
1.0,
0.8,
0.6,
0.4,
0.2,
0.08
]

def Spec2RGB(SpecType):
  if not _StarSpec_.has_key(SpecType):
    SpecType='-'
  return _StarSpec_[SpecType]

def Vmag2Size(Vmag):
  ret=_StarSize_[5]
  if Vmag<1: ret=_StarSize_[0]
  else:
    if Vmag<2: ret=_StarSize_[1]
    else:
      if Vmag<3: ret=_StarSize_[2]
      else:
        if Vmag<4: ret=_StarSize_[3]
        else:
         if Vmag<5: ret=_StarSize_[4]
  return ret