Compute local solar time for given (lat, lon)
338153
def solarTime(utc_dt, lat, lon):
    
    import ephem
    o = ephem.Observer()
    o.date = utc_dt
    o.lat = str(lat)
    o.lon = str(lon)
    sun = ephem.Sun()
    sun.compute(o)
    hour_angle = o.sidereal_time() - sun.ra
    rad = str(ephem.hours(hour_angle + ephem.hours('12:00')).norm)
    t = datetime.strptime(rad, '%H:%M:%S.%f')
    solar_dt = datetime.combine(utc_dt.date(), t.time()) 
    return solar_dt
1126277
def getTimeZone(lat, lon):
    
    #Need to fix for Python 2.x and 3.X support
    import urllib.request, urllib.error, urllib.parse
    import xml.etree.ElementTree as ET
    #http://api.askgeo.com/v1/918/aa8292ec06199d1207ccc15be3180213c984832707f0cbf3d3859db279b4b324/query.xml?points=37.78%2C-122.42%3B40.71%2C-74.01&databases=Point%2CTimeZone%2CAstronomy%2CNaturalEarthCountry%2CUsState2010%2CUsCounty2010%2CUsCountySubdivision2010%2CUsTract2010%2CUsBlockGroup2010%2CUsPlace2010%2CUsZcta2010
    req = "http://api.askgeo.com/v1/918/aa8292ec06199d1207ccc15be3180213c984832707f0cbf3d3859db279b4b324/query.xml?points="+str(lat)+"%2C"+str(lon)+"&databases=TimeZone"
    opener = urllib.request.build_opener()
    f = opener.open(req)
    tree = ET.parse(f)
    root = tree.getroot()
    #Check response
    tzid = None
    if root.attrib['code'] == '0':
        tz = list(root.iter('TimeZone'))[0]
        #shortname = tz.attrib['ShortName']
        tzid = tz.attrib['TimeZoneId']
    return tzid
1193786
def coord_from_area(self, x, y, lat, lon, width, ground_width):
		''''''

		pixel_width = ground_width / float(width)
		dx = x * pixel_width
		dy = y * pixel_width

		return mp_util.gps_offset(lat, lon, dx, -dy)
69623
def G(self, ID, lat, lon):
        
        
        # Equatorial coordinates
        eqM = utils.eqCoords(lon, lat)
        eqZ = eqM
        if lat != 0:
            eqZ = utils.eqCoords(lon, 0)
        
        return {
            'id': ID,
            'lat': lat,
            'lon': lon,
            'ra': eqM[0],
            'decl': eqM[1],
            'raZ': eqZ[0],
            'declZ': eqZ[1],
        }
712087
def nextSunset(jd, lat, lon):
    
    return swe.sweNextTransit(const.SUN, jd, lat, lon, 'SET')
465654
function phase (phase_date) {
  if (!phase_date) {
    phase_date = new Date()
  }
  phase_date = julian.fromDate(phase_date)

  const day = phase_date - EPOCH

  // calculate sun position
  const sun_mean_anomaly =
    (360.0 / 365.2422) * day +
    (ECLIPTIC_LONGITUDE_EPOCH - ECLIPTIC_LONGITUDE_PERIGEE)
  const sun_true_anomaly =
    2 * todeg(Math.atan(
      Math.sqrt((1.0 + ECCENTRICITY) / (1.0 - ECCENTRICITY)) *
      Math.tan(0.5 * kepler(sun_mean_anomaly, ECCENTRICITY))
    ))
  const sun_ecliptic_longitude =
    ECLIPTIC_LONGITUDE_PERIGEE + sun_true_anomaly
  const sun_orbital_distance_factor =
    (1 + ECCENTRICITY * dcos(sun_true_anomaly)) /
    (1 - ECCENTRICITY * ECCENTRICITY)

  // calculate moon position
  const moon_mean_longitude =
    MOON_MEAN_LONGITUDE_EPOCH + 13.1763966 * day
  const moon_mean_anomaly =
    moon_mean_longitude - 0.1114041 * day - MOON_MEAN_PERIGEE_EPOCH
  const moon_evection =
    1.2739 * dsin(
      2 * (moon_mean_longitude - sun_ecliptic_longitude) - moon_mean_anomaly
    )
  const moon_annual_equation =
    0.1858 * dsin(sun_mean_anomaly)
  // XXX: what is the proper name for this value?
  const moon_mp =
    moon_mean_anomaly +
    moon_evection -
    moon_annual_equation -
    0.37 * dsin(sun_mean_anomaly)
  const moon_equation_center_correction =
    6.2886 * dsin(moon_mp)
  const moon_corrected_longitude =
    moon_mean_longitude +
    moon_evection +
    moon_equation_center_correction -
    moon_annual_equation +
    0.214 * dsin(2.0 * moon_mp)
  const moon_age =
    fixangle(
      moon_corrected_longitude -
      sun_ecliptic_longitude +
      0.6583 * dsin(
        2 * (moon_corrected_longitude - sun_ecliptic_longitude)
      )
    )
  const moon_distance =
    (MOON_SMAXIS * (1.0 - MOON_ECCENTRICITY * MOON_ECCENTRICITY)) /
    (1.0 + MOON_ECCENTRICITY * dcos(moon_mp + moon_equation_center_correction))

  return {
    phase: (1.0 / 360.0) * moon_age,
    illuminated: 0.5 * (1.0 - dcos(moon_age)),
    age: (SYNODIC_MONTH / 360.0) * moon_age,
    distance: moon_distance,
    angular_diameter: MOON_ANGULAR_SIZE_SMAXIS / moon_distance,
    sun_distance: SUN_SMAXIS / sun_orbital_distance_factor,
    sun_angular_diameter: SUN_ANGULAR_SIZE_SMAXIS * sun_orbital_distance_factor
  }
}
642032
def earth_sun_d(dtime):
        
        doy = int(dtime.strftime('%j'))
        rad_term = 0.9856 * (doy - 4) * pi / 180
        distance_au = 1 - 0.01672 * cos(rad_term)
        return distance_au
584920
def datetimes(self):
        """"""
        if self.timestep == 1:
            return tuple(dt.add_minute(30) for dt in
                         self.direct_normal_irradiance.datetimes)
        else:
            return self.direct_normal_irradiance.datetimes
100470
def get_alt_az(utc_time, lon, lat):
    
    lon = np.deg2rad(lon)
    lat = np.deg2rad(lat)

    ra_, dec = sun_ra_dec(utc_time)
    h__ = _local_hour_angle(utc_time, lon, ra_)
    return (np.arcsin(np.sin(lat) * np.sin(dec) +
                      np.cos(lat) * np.cos(dec) * np.cos(h__)),
            np.arctan2(-np.sin(h__), (np.cos(lat) * np.tan(dec) -
                                      np.sin(lat) * np.cos(h__))))
72612
def _analemma_suns(self):
        
        for h in xrange(0, 24):
            if self._analemma_position(h) < 0:
                continue
            elif self._analemma_position(h) == 0:
                chours = []
                # this is an hour that not all the hours are day or night
                prevhour = self.latitude <= 0
                num_of_days = 8760 if not self.is_leap_year else 8760 + 24
                for hoy in xrange(h, num_of_days, 24):
                    thishour = self.calculate_sun_from_hoy(hoy).is_during_day
                    if thishour != prevhour:
                        if not thishour:
                            hoy -= 24
                        dt = DateTime.from_hoy(hoy, self.is_leap_year)
                        chours.append((dt.month, dt.day, dt.hour))
                    prevhour = thishour
                tt = []
                for hcount in range(int(len(chours) / 2)):
                    st = chours[2 * hcount]
                    en = chours[2 * hcount + 1]
                    if self.latitude >= 0:
                        tt = [self.calculate_sun(*st)] + \
                            [self.calculate_sun(st[0], d, h)
                             for d in xrange(st[1] + 1, 29, 7)] + \
                            [self.calculate_sun(m, d, h)
                             for m in xrange(st[0] + 1, en[0])
                             for d in xrange(3, 29, 7)] + \
                            [self.calculate_sun(en[0], d, h)
                             for d in xrange(3, en[1], 7)] + \
                            [self.calculate_sun(*en)]
                    else:
                        tt = [self.calculate_sun(*en)] + \
                            [self.calculate_sun(en[0], d, h)
                             for d in xrange(en[1] + 1, 29, 7)] + \
                            [self.calculate_sun(m, d, h) for m in xrange(en[0] +
                                                                         1, 13)
                             for d in xrange(3, 29, 7)] + \
                            [self.calculate_sun(m, d, h) for m in xrange(1, st[0])
                             for d in xrange(3, 29, 7)] + \
                            [self.calculate_sun(st[0], d, h)
                             for d in xrange(3, st[1], 7)] + \
                            [self.calculate_sun(*st)]
                    yield tt
            else:
                yield tuple(self.calculate_sun((m % 12) + 1, d, h)
                            for m in xrange(0, 13) for d in (7, 14, 21))[:-2]
725224
def mktime_tz(data):
    """"""
    if data[9] is None:
        # No zone info, so localtime is better assumption than GMT
        return time.mktime(data[:8] + (-1,))
    else:
        t = calendar.timegm(data)
        return t - data[9]
300362
def nextSolarReturn(date, lon):
    
    jd = eph.nextSolarReturn(date.jd, lon)
    return Datetime.fromJD(jd, date.utcoffset)
207176
public Moonlight on(CalendarDate date) {

        // initialization
        PlainDate d = SolarTime.toGregorian(date);
        Timezone tz = Timezone.of(this.observerZoneID);
        Moment start =
            ((tz.getHistory() == null)
                ? d.at(PlainTime.midnightAtStartOfDay()).in(tz)
                : d.atFirstMoment(this.observerZoneID));

        double mjd0 = JulianDay.ofMeanSolarTime(start).getMJD();
        double longitudeRad = Math.toRadians(this.longitude);
        double cosLatitude = Math.cos(Math.toRadians(this.latitude));
        double sinLatitude = Math.sin(Math.toRadians(this.latitude));
        double geodeticAngle = StdSolarCalculator.TIME4J.getGeodeticAngle(this.latitude, this.altitude);
        double refraction = AstroUtils.refractionFactorOfStdAtmosphere(this.altitude) * SolarTime.STD_REFRACTION / 60;
        double deltaT = TimeScale.deltaT(d);
        double hour = 1.0;
        double y_minus = sinAlt(mjd0, 0.0, longitudeRad, cosLatitude, sinLatitude, geodeticAngle, refraction, deltaT);
        double[] result = new double[4];

        // declaration of result data
        boolean above = (y_minus > 0.0); // at start of day
        boolean rises = false;
        boolean sets = false;
        double risingHour = Double.NaN;
        double settingHour = Double.NaN;

        // loop over 2-hour-search-intervals applying quadratic interpolation
        do {
            double y_0 =
                sinAlt(mjd0, hour, longitudeRad, cosLatitude, sinLatitude, geodeticAngle, refraction, deltaT);
            double y_plus =
                sinAlt(mjd0, hour + 1, longitudeRad, cosLatitude, sinLatitude, geodeticAngle, refraction, deltaT);
            int count =
                interpolate(y_minus, y_0, y_plus, result);
            if (count == 1) {
                double root = result[2];
                if (Double.isNaN(root)) {
                    root = result[3];
                }
                if (y_minus < 0.0) {
                    risingHour = hour + root;
                    rises = true;
                } else {
                    settingHour = hour + root;
                    sets = true;
                }
            } else if (count == 2) {
                if (result[1] < 0.0) {
                    risingHour = hour + result[3];
                    settingHour = hour + result[2];
                } else {
                    risingHour = hour + result[2];
                    settingHour = hour + result[3];
                }
                rises = true;
                sets = true;
            }
            y_minus = y_plus;
            hour += 2.0;
        } while (!((hour > 25.0) || (rises && sets))); // (> 25.0)-condition cares about possible 25-h-day (end-of-DST)

        // evaluate moonrise and moonset
        Moment rising = null;
        Moment setting = null;

        if (rises) {
            rising = add(start, risingHour);
            if (!rising.toZonalTimestamp(this.observerZoneID).getCalendarDate().equals(d)) {
                rising = null;
                // rises = false;
            }
        }
        if (sets) {
            setting = add(start, settingHour);
            if (!setting.toZonalTimestamp(this.observerZoneID).getCalendarDate().equals(d)) {
                setting = null;
                // sets = false;
            }
        }

        return new Moonlight(d, this.observerZoneID, rising, setting, above);

    }
987029
def observer_position(time, lon, lat, alt):
    

    lon = np.deg2rad(lon)
    lat = np.deg2rad(lat)

    theta = (gmst(time) + lon) % (2 * np.pi)
    c = 1 / np.sqrt(1 + F * (F - 2) * np.sin(lat)**2)
    sq = c * (1 - F)**2

    achcp = (A * c + alt) * np.cos(lat)
    x = achcp * np.cos(theta)  # kilometers
    y = achcp * np.sin(theta)
    z = (A * sq + alt) * np.sin(lat)

    vx = -MFACTOR * y  # kilometers/second
    vy = MFACTOR * x
    vz = 0

    return (x, y, z), (vx, vy, vz)
35836
def eqCoords(lon, lat):
    
    # Convert to radians
    _lambda = math.radians(lon)
    _beta = math.radians(lat)
    _epson = math.radians(23.44)  # The earth's inclination
    
    # Declination in radians
    decl = math.asin(math.sin(_epson) * math.sin(_lambda) * math.cos(_beta) + \
           math.cos(_epson) * math.sin(_beta))
    
    # Equatorial Distance in radians
    ED = math.acos(math.cos(_lambda) * math.cos(_beta) / math.cos(decl))
    
    # RA in radians
    ra = ED if lon < 180 else math.radians(360) - ED
    
    # Correctness of RA if longitude is close to 0º or 180º in a radius of 5º
    if (abs(angle.closestdistance(lon, 0)) < 5 or 
        abs(angle.closestdistance(lon, 180)) < 5):
            a = math.sin(ra) * math.cos(decl)
            b = math.cos(_epson) * math.sin(_lambda) * math.cos(_beta) - \
                math.sin(_epson) * math.sin(_beta)
            if (math.fabs(a-b) > 0.0003):
                ra = math.radians(360) - ra
                
    return (math.degrees(ra), math.degrees(decl))
430330
def solarzenithangle(time: datetime, glat: float, glon: float, alt_m: float) -> tuple:
    
    time = totime(time)

    obs = EarthLocation(lat=glat*u.deg, lon=glon*u.deg, height=alt_m*u.m)
    times = Time(time, scale='ut1')
    sun = get_sun(times)
    sunobs = sun.transform_to(AltAz(obstime=times, location=obs))

    return 90 - sunobs.alt.degree, sun, sunobs
606553
def _calculate_sunrise_hour_angle(self, solar_dec, depression=0.833):
        """"""

        hour_angle_arg = math.degrees(math.acos(
            math.cos(math.radians(90 + depression)) /
            (math.cos(math.radians(self.latitude)) * math.cos(
                math.radians(solar_dec))) -
            math.tan(math.radians(self.latitude)) *
            math.tan(math.radians(solar_dec))
        ))

        return hour_angle_arg