Fluxion with the potential energy of the earth-sun sytem
973910
def make_force_ES(q_vars, mass):
    """"""
    # Build the potential energy fluxion; just one pair of bodies
    U =  U_ij(q_vars, mass, 0, 1)
    # Varname arrays for both the coordinate system and U
    vn_q = np.array([q.var_name for q in q_vars])
    vn_fl = np.array(sorted(U.var_names))
    # Permutation array for putting variables in q in the order expected by U (alphabetical)
    q2fl = np.array([np.argmax((vn_q == v)) for v in vn_fl])
    # Permutation array for putting results of U.diff() in order of q_vars
    fl2q = np.array([np.argmax((vn_fl == v)) for v in vn_q])
    # Return a force function from this potential
    force_func = lambda q: -U.diff(q[q2fl]).squeeze()[fl2q]
    return force_func
157731
def expected_mag(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    m_body = Vector3(RAW_IMU.xmag, RAW_IMU.ymag, RAW_IMU.zmag)
    field_strength = m_body.length()

    m = rotation(ATTITUDE)

    r = Matrix3()
    r.from_euler(0, -radians(inclination), radians(declination))
    m_earth = r * Vector3(field_strength, 0, 0)

    return m.transposed() * m_earth
1146758
def get_potential_energy(self, a):
        """"""
        e = 0.0
        for c in self.calcs:
            e += c.get_potential_energy(a)
        return e
704600
def _compute_heating_rates(self):
        ''''''
        self.LW_to_atm = self._compute_emission()
        self.SW_to_atm = self._compute_reflected_flux()
        self.heating_rate['Ts'] = ( self.LW_from_atm - self.LW_to_atm
                                  + self.SW_from_atm - self.SW_to_atm )
1354590
def _get_flux_bounds(self, r_id, model, flux_limits, equation):
        
        if r_id not in flux_limits or flux_limits[r_id][0] is None:
            if equation.direction == Direction.Forward:
                lower = 0
            else:
                lower = -model.default_flux_limit
        else:
            lower = flux_limits[r_id][0]

        if r_id not in flux_limits or flux_limits[r_id][1] is None:
            if equation.direction == Direction.Reverse:
                upper = 0
            else:
                upper = model.default_flux_limit
        else:
            upper = flux_limits[r_id][1]

        if lower % 1 == 0:
            lower = int(lower)
        if upper % 1 == 0:
            upper = int(upper)
        return text_type(lower), text_type(upper)
370173
function forward$2(p) {
        var lon = p.x;
        var lat = p.y;

        var delta_lon = adjust_lon(lon - this.long0);
        var con;
        var x, y;
        var sin_phi = Math.sin(lat);
        var cos_phi = Math.cos(lat);

        if (!this.es) {
            var b = cos_phi * Math.sin(delta_lon);

            if ((Math.abs(Math.abs(b) - 1)) < EPSLN) {
                return (93);
            }
            else {
                x = 0.5 * this.a * this.k0 * Math.log((1 + b) / (1 - b)) + this.x0;
                y = cos_phi * Math.cos(delta_lon) / Math.sqrt(1 - Math.pow(b, 2));
                b = Math.abs(y);

                if (b >= 1) {
                    if ((b - 1) > EPSLN) {
                        return (93);
                    }
                    else {
                        y = 0;
                    }
                }
                else {
                    y = Math.acos(y);
                }

                if (lat < 0) {
                    y = -y;
                }

                y = this.a * this.k0 * (y - this.lat0) + this.y0;
            }
        }
        else {
            var al = cos_phi * delta_lon;
            var als = Math.pow(al, 2);
            var c = this.ep2 * Math.pow(cos_phi, 2);
            var cs = Math.pow(c, 2);
            var tq = Math.abs(cos_phi) > EPSLN ? Math.tan(lat) : 0;
            var t = Math.pow(tq, 2);
            var ts = Math.pow(t, 2);
            con = 1 - this.es * Math.pow(sin_phi, 2);
            al = al / Math.sqrt(con);
            var ml = pj_mlfn(lat, sin_phi, cos_phi, this.en);

            x = this.a * (this.k0 * al * (1 +
                als / 6 * (1 - t + c +
                als / 20 * (5 - 18 * t + ts + 14 * c - 58 * t * c +
                als / 42 * (61 + 179 * ts - ts * t - 479 * t))))) +
                this.x0;

            y = this.a * (this.k0 * (ml - this.ml0 +
                sin_phi * delta_lon * al / 2 * (1 +
                als / 12 * (5 - t + 9 * c + 4 * cs +
                als / 30 * (61 + ts - 58 * t + 270 * c - 330 * t * c +
                als / 56 * (1385 + 543 * ts - ts * t - 3111 * t)))))) +
                this.y0;
        }

        p.x = x;
        p.y = y;

        return p;
    }
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
1319644
def solar_constant(self):
        """"""
        if self.wavenumber is not None:
            return np.trapz(self.irradiance, self.wavenumber)
        elif self.wavelength is not None:
            return np.trapz(self.irradiance, self.wavelength)
        else:
            raise TypeError('Neither wavelengths nor wavenumbers available!')
794900
def earth_accel2_df(IMU,IMU2,ATT):
    ''''''
    r = rotation_df(ATT)
    accel1 = Vector3(IMU.AccX, IMU.AccY, IMU.AccZ)
    accel2 = Vector3(IMU2.AccX, IMU2.AccY, IMU2.AccZ)
    accel = 0.5 * (accel1 + accel2)
    return r * accel
163781
public function getDistance($lat1, $lgn1, $lat2, $lgn2){
        $earth_radius = 6371;
        $alpha = $lat1 - $lat2;
        $beta = $lgn1 - $lgn2;
        $h = pow((sin(deg2rad($alpha/2))), 2) +  cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * pow(sin(deg2rad($beta/2)),2);
        $dist = 2*$earth_radius*sqrt($h);

        return $dist;
    }
108560
def configuration_ES(t0: date, t1: Optional[date] = None, 
                     steps_per_day: int = None) -> Tuple[np.ndarray, np.ndarray]:
    
    # Default steps_per_day = 1
    if steps_per_day is None:
        steps_per_day = 1
        
    # Time step dt is 1.0 over steps per day
    dt: float = 1.0 / float(steps_per_day)

    # Default t1 to one day after t0
    if t1 is not None:
        # Convert t to a julian day
        jd0: int = julian_day(t0)
        jd1: int = julian_day(t1)
    else:
        jd0: int = julian_day(t0)
        jd1: int = jd0 + dt
    # Pass the times as an array of julian days
    jd: np.ndarray = np.arange(jd0, jd1, dt)

    # Position and velocity of the sun as arrays of length 3
    sun_id: int = jpl_body_id['sun']
    pos_sun, vel_sun = jpl_kernel[0, sun_id].compute_and_differentiate(jd)

    # Position and velocity of the earth as arrays of length 3
    earth_id: int = jpl_body_id['earth']
    pos_earth, vel_earth = jpl_kernel[0, earth_id].compute_and_differentiate(jd)

    # Convert positions from km to meters (multiply by km2m)
    q = np.vstack([pos_sun, pos_earth]).T * km2m
    # Convert velocities from km / day to meters / sec (multiply by km2m, divide by day2sec)
    v = np.vstack([vel_sun, vel_earth]).T * (km2m / day2sec)

    # Return tuple of Tx6 arrays for position q and velocity v
    return q, v
845145
def _total_flux(self, r_eff, I_eff, n_sersic):
        
        bn = self.b_n(n_sersic)
        return I_eff * r_eff**2 * 2 * np.pi * n_sersic * np.exp(bn) / bn**(2*n_sersic) * scipy.special.gamma(2*n_sersic)
964984
def major_flux(self, fraction=0.9):
        

        """
        (paths, pathfluxes) = self.pathways(fraction=fraction)
        return self._pathways_to_flux(paths, pathfluxes, n=self.nstates)
1001309
def _gti_dirint_gte_90(poa_global, aoi, solar_zenith, solar_azimuth,
                       surface_tilt, times, kt_prime,
                       pressure=101325., temp_dew=None, albedo=.25):
    
    kt_prime_gte_90 = _gti_dirint_gte_90_kt_prime(aoi, solar_zenith,
                                                  solar_azimuth, times,
                                                  kt_prime)

    I0 = get_extra_radiation(times, 1370, 'spencer')
    airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
    airmass = atmosphere.get_absolute_airmass(airmass, pressure)
    kt = kt_prime_gte_90 * _kt_kt_prime_factor(airmass)
    disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)

    dni_gte_90 = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
                                          False, temp_dew)

    dni_gte_90_proj = dni_gte_90 * tools.cosd(solar_zenith)
    cos_surface_tilt = tools.cosd(surface_tilt)

    # isotropic sky plus ground diffuse
    dhi_gte_90 = (
        (2 * poa_global - dni_gte_90_proj * albedo * (1 - cos_surface_tilt)) /
        (1 + cos_surface_tilt + albedo * (1 - cos_surface_tilt)))

    ghi_gte_90 = dni_gte_90_proj + dhi_gte_90

    return ghi_gte_90, dni_gte_90, dhi_gte_90
821214
def evaluate(self, x, *args):
        
        a = (self.amplitude * np.ones_like(x)) * self._flux_unit
        y = units.convert_flux(x, a, units.PHOTLAM)
        return y.value
972523
def sidereal_time(t):
    """"""

    # Compute the Earth Rotation Angle.  Time argument is UT1.

    theta = earth_rotation_angle(t.ut1)

    # The equinox method.  See Circular 179, Section 2.6.2.
    # Precession-in-RA terms in mean sidereal time taken from third
    # reference, eq. (42), with coefficients in arcseconds.

    t = (t.tdb - T0) / 36525.0
    st =        ( 0.014506 +
        (((( -    0.0000000368   * t
             -    0.000029956  ) * t
             -    0.00000044   ) * t
             +    1.3915817    ) * t
             + 4612.156534     ) * t)

    # Form the Greenwich sidereal time.

    return (st / 54000.0 + theta * 24.0) % 24.0
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