text:
# -*- coding: utf-8 -*-

import pytest

import numpy as np

from threedidepth import morton

def test_morton_array():
    result = morton.morton_array((8, 8)).tolist()
    expected = [[ 0,  1,  4,  5, 16, 17, 20, 21],  # noqa: E201
                [ 2,  3,  6,  7, 18, 19, 22, 23],  # noqa: E201
                [ 8,  9, 12, 13, 24, 25, 28, 29],  # noqa: E201
                [10, 11, 14, 15, 26, 27, 30, 31],
                [32, 33, 36, 37, 48, 49, 52, 53],
                [34, 35, 38, 39, 50, 51, 54, 55],
                [40, 41, 44, 45, 56, 57, 60, 61],
                [42, 43, 46, 47, 58, 59, 62, 63]]
    assert result == expected
    assert morton.morton_array((2,)).tolist() == [0, 1]
    assert morton.morton_array((1, 2)).tolist() == [[0, 1]]
    with pytest.raises(ValueError):
        morton.morton_array((2 ** 32, 2 ** 32))

def test_group():
    array = np.array([0.6, 0.6, 0.6, 0.5, 0.7, 0.5, 0.7, 0.5, 0.7])
    expected = [[3, 5, 7], [0, 1, 2], [4, 6, 8]]
    result = [i.tolist() for i in morton.group(array)]
    assert result == expected

def test_rasterize():
    # single point
    points = np.array([[0, 0]])
    expected = [[0]]
    result = morton.rasterize(points)[0].tolist()
    assert result == expected

    # horizontal
    points = np.array([[0, 0], [1, 0]])
    expected = [[0, 1]]
    result = morton.rasterize(points)[0].tolist()
    assert result == expected

    # vertical
    points = np.array([[0, 0], [1, 0]])
    expected = [[0, 1]]
    result = morton.rasterize(points)[0].tolist()
    assert result == expected

    # quadtree
    points = np.array([[1, 1], [3, 1], [1, 3], [3, 3], [6, 2]])
    expected = [
        [2, 3, 4],
        [0, 1, 5],
    ]
    result = morton.rasterize(points)[0].tolist()
    assert result == expected

def test_reorder():
    points = np.array([[1, 3], [3, 1], [6, 2], [3, 3], [1, 1]])
    s1 = np.array([1, 2, 3, 4, 5])

    result_points, result_s1 = morton.reorder(points, s1)

    expected_points = [[1, 3], [3, 3], [1, 1], [3, 1], [6, 2]]
    expected_s1 = [1, 4, 5, 2, 3]

    assert result_points.tolist() == expected_points
    assert result_s1.tolist() == expected_s1

max_stars_repo_name: nens/threedidepth

max_issues_repo_name: nens/threedidepth

max_forks_repo_name: nens/threedidepth


text:
from sklearn.decomposition import NMF
import numpy as np
import pandas as pd

# 1. Matrix Factorization: NMF
print('========= 1. Matrix Factorization using NMF =========')
ratings = pd.read_csv('chapter04/data/ml-latest-small/ratings.csv')
ratings = pd.pivot_table(
    ratings,
    index='userId',
    columns='movieId'
)
ratings = ratings.fillna(0)
print(ratings.shape)

nmf = NMF(n_components=10)
p = nmf.fit_transform(ratings)
q = nmf.components_
r_hat = np.dot(p, q)

print(p.shape)
print(q.shape)
print(r_hat.shape)

# 2. Matrix Factorization: SVD
print('========= 2. Matrix Factorization using SVD =========')
sample_matrix = np.random.uniform(0, 5, size=(9, 6))
print('Original Matrix')
print(sample_matrix)

U, S, V = np.linalg.svd(sample_matrix)
print(U.shape)
print(S.shape)
print(V.shape)

k = 2
print(U[:, :k] @ np.diag(S[:k]) @ V[:k, :])

ext: py

reason: import numpy


text:
import collections
import numpy as np
from scipy import spatial as sp
from .base import *
from .data import *
from .bse import *

class Core(DictNpArrayMix):
    """ Core of a star cluster
    Keys: (class members)
        time (1D): evolved time of the star cluster
        pos  (2D,3): center of the star cluster
        vel  (2D,3): center velocity of the cluster
        rc   (1D): core radius of the star cluster
    """

    def __init__(self, _dat=None, _offset=int(0), _append=False, **kwargs):
        """ DictNpArrayMix type initialzation, see help(DictNpArrayMix.__init__)
        """
        keys  = [['time',np.float64],['pos', (np.float64,3)],['vel', (np.float64,3)], ['rc', np.float64]]
        DictNpArrayMix.__init__(self, keys, _dat, _offset, _append, **kwargs)

    def calcPotentialCenter(self, single, binary, G):
        """ Calculate potential center of the system and return the result
        r_cm = \sum_i pot_i *r_i /\sum_i pot_i (only count pot_i <0)

        Parameters
        ----------
        single: inhermited SimpleParticle
            Single particle data
        binary: Binary
            Binary particle data
        G: float
            Gravitational constant

        Return
        ----------
        cm_pos: c.m. position, numpy.ndarray with shape (*,3)
        cm_vel: c.m. velocity, numpy.ndarray with shape (*,3)

        """
        pot_s_sel = (single.pot<0)
        pot_s = single.pot[pot_s_sel]
        pos_s = single.pos[pot_s_sel]
        vel_s = single.vel[pot_s_sel]

        binary.calcPot(G)
        pot_b = binary.pot
        pos_b = binary.pos
        vel_b = binary.vel

        pot_sum = pot_s.sum() + pot_b.sum()
        cm_pos = np.array([(np.sum(pot_s*pos_s[:,i])+np.sum(pot_b*pos_b[:,i]))/pot_sum for i in range(3)])
        cm_vel = np.array([(np.sum(pot_s*vel_s[:,i])+np.sum(pot_b*vel_b[:,i]))/pot_sum for i in range(3)])

        self.pos = np.append(self.pos,[cm_pos],axis=0)
        self.vel = np.append(self.vel,[cm_vel],axis=0)

        return cm_pos, cm_vel

    def calcDensityAndCenter(self, particle, kdtree):
        """ Calculate density based on nearest six neighbors and return the result (Casertano & Hut 1985)

        Parameters
        ----------
        particle: inherited SimpleParticle
            Particle data set
        kdtree: scipy.spatial.cKDtree
            3D KDTree of all particle positions, can be generated from findPair

        Return
        ----------
        cm_pos: c.m. position, numpy.ndarray with shape (*,3)
        cm_vel: c.m. velocity, numpy.ndarray with shape (*,3)
        """
        # 6 nearest neighbors
        nb_r_list6, nb_index_list6 = kdtree.query(particle.pos,k=6)
        nb_mass_tot6=np.sum(particle[nb_index_list6].mass,axis=1) + particle.mass

        nb_inv_r6 = 1/nb_r_list6[:,5]
        rho = nb_mass_tot6*(nb_inv_r6*nb_inv_r6*nb_inv_r6)
        particle.addNewMember('density',rho)
        rho_tot = rho.sum()

        cm_pos = np.array([np.sum(rho*particle.pos[:,i])/rho_tot for i in range(3)])
        cm_vel = np.array([np.sum(rho*particle.vel[:,i])/rho_tot for i in range(3)])
        self.pos = np.append(self.pos, [cm_pos], axis=0)
        self.vel = np.append(self.vel, [cm_vel], axis=0)

        return cm_pos, cm_vel

    def calcCoreRadius(self, particle):
        """ Calculate core radius, using Casertano & Hut (1985) method:
        rc = sqrt(\sum_i rho_i^2 r_i^2 / \sum_i rho_i^2)

        Parameters
        ----------
        particle: inherited SimpleParticle
            Particle data set

        Return
        ----------
        rc: core radius

        """
        rho2 = particle.density*particle.density
        rc = np.sqrt((particle.r2*rho2).sum()/(rho2.sum()))
        self.rc = np.append(self.rc, rc)

        return rc

    def addTime(self, time):
        """ Append a new time to current member 'time'
        """
        self.time = np.append(self.time, time)

    def toSkyCoord(self, **kwargs):
        """ generate astropy.coordinates.SkyCoord data in galactocentric frame

        Parameters
        -----------------
        kwargs: dict()
            pos_unit: astropy.units (units.pc)
                 position unit of the particle data
            vel_unit: astropy.units (units.pc/units.Myr)
                 velocity unit of the particle data
            galcen_distance: floating with length units (8.0*units.kpc [Galpy])
                 galactic central distance of the Sun
            z_sun: floating with length units (15.0*units.pc [Galpy])
                 z direction distance of the Sun
            galcen_v_sun: astropy.coordinates.CartesianDifferential ([10.0, 235.0, 7.0]*units.km/units.s [Galpy])
                 velocity of the Sun

        Return
        ----------------
        core_g: astropy.coordinates.SkyCoord
            core c.m. data using SkyCoord
        """
        import astropy
        from astropy.coordinates import SkyCoord  # High-level coordinates
        from astropy.coordinates import ICRS, Galactic, Galactocentric, FK4, FK5  # Low-level frames
        from astropy.coordinates import Angle, Latitude, Longitude  # Angles
        from astropy.coordinates import CartesianDifferential
        import astropy.units as u

        pos_unit = u.pc
        if ('pos_unit' in kwargs.keys()): pos_unit = kwargs['pos_unit']
        vel_unit = u.pc/u.Myr
        if ('vel_unit' in kwargs.keys()): vel_unit = kwargs['vel_unit']

        parameters={'galcen_distance':8.0*u.kpc, 'z_sun':15.*u.pc, 'galcen_v_sun':CartesianDifferential([10.0,235.,7.]*u.km/u.s)}
        for key in parameters.keys():
            if key in kwargs.keys():
                parameter[key] = kwargs[key]

        core_g = SkyCoord(x=self.pos[:,0]*pos_unit,
                          y=self.pos[:,1]*pos_unit,
                          z=self.pos[:,2]*pos_unit,
                          v_x=self.vel[:,0]*vel_unit,
                          v_y=self.vel[:,1]*vel_unit,
                          v_z=self.vel[:,2]*vel_unit,
                          frame='galactocentric', representation_type='cartesian', **parameters)
        return core_g

class LagrangianVelocity(DictNpArrayMix):
    """ Lagrangian velocity component
        Each velocity member is a 2D numpy.ndarray, the row contain the values corresponding to each Langragian radius and the core radius (last value).
        For example, in default case, mass_fraction array is np.array([0.1, 0.3, 0.5, 0.7, 0.9]).
        The member 'abs' has the shape (*,6), each row contain the velocity referring to the 0.1, 0.3, 0.5, 0.7, 0.9 of Langragian radii and the core radius (last value).

    Keys: (class members)
        abs (2D,n_frac): absolute 3D velocity value
        x   (2D,n_frac): velocity in x axis
        y   (2D,n_frac): velocity in y axis
        z   (2D,n_frac): velocity in z axis
        rad (2D,n_frac): velocity in radial direction
        tan (2D,n_frac): velocity in tangential direction
        rot (2D,n_frac): rotational velocity in x-y plane

        n_frac: number of Lagrangian radii + 1, determined by keyword argument 'mass_fraction'
    """

    def __init__(self, _dat=None, _offset=int(0), _append=False, **kwargs):
        """ DictNpArrayMix type initialzation, see help(DictNpArrayMix.__init__)

        Parameters
        ----------
        keyword arguments:
            mass_fraction: 1D numpy.ndarray (np.array([0.1, 0.3, 0.5, 0.7, 0.9]))
                An array to indicate the mass fractions to calculate lagrangian radii.
        """

        m_frac=np.array([0.1,0.3,0.5,0.7,0.9])
        if ('mass_fraction' in kwargs.keys()): m_frac=kwargs['mass_fraction'].copy()
        n_frac = m_frac.size + 1
        keys  = [['abs',(np.float64,n_frac)],['x',(np.float64,n_frac)],['y',(np.float64,n_frac)],['z',(np.float64,n_frac)],['rad',(np.float64,n_frac)],['tan',(np.float64,n_frac)],['rot',(np.float64,n_frac)]] # all, x, y, z, radial, tangential, rotational
        DictNpArrayMix.__init__(self, keys, _dat, _offset, _append, **kwargs)
        self.initargs['mass_fraction'] = m_frac

class Lagrangian(DictNpArrayMix):
    """ Lagrangian parameters
    Keys: (class members)
        r (2D,n_frac): Lagrangian radii
        m (2D,n_frac): average mass referring to Lagrangian radii
        n (2D,n_frac): number of particles referring to Lagrangian radii
        vel   (LagrangianVelocity): average velocity referring to Lagrangian radii
        sigma (LagrangianVelocity): velocity dispersion referring to Lagrangian radii

        Additional members ff the keyword argument 'calc_energy' is true:
          epot (2D, n_frac): potential energy
          vr (2D, n_frac); virial ratio -(2*ekin/epot); notice that here the external potential energy is reduced by half.
          *epot_ext (2D, n_frac): external potential energy (exists if the keyword argument 'external_mode' is switch on)
    """

    def __init__(self, _dat=None, _offset=int(0), _append=False, **kwargs):
        """ DictNpArrayMix type initialzation, see help(DictNpArrayMix.__init__)

        Parameters
        ----------
        keyword arguments:
            mass_fraction: 1D numpy.ndarray (np.array([0.1, 0.3, 0.5, 0.7, 0.9]))
                An array to indicate the mass fractions to calculate lagrangian radii.
            calc_energy: bool (False)
                If true, add kinetic, potential energies and virial ratio calculation
            external_mode: string (none)
                PeTar external mode (set in configure): galpy, none
                If it is not none, epot_ext will be added if calc_energy=True
            calc_multi_rc: bool (False)
                If true, calculate core radius using KDTree and correct the central position
        """
        m_frac=np.array([0.1,0.3,0.5,0.7,0.9])
        if ('mass_fraction' in kwargs.keys()): m_frac=kwargs['mass_fraction'].copy()
        calc_energy=False
        if ('calc_energy' in kwargs.keys()): calc_energy=kwargs['calc_energy']
        external_mode='none'
        if ('external_mode' in kwargs.keys()): external_mode=kwargs['external_mode']
        calc_multi_rc=False
        if ('calc_multi_rc' in kwargs.keys()): calc_multi_rc=kwargs['calc_multi_rc']

        n_frac = m_frac.size + 1
        keys = [['r', (np.float64, n_frac)],['m', (np.float64, n_frac)],['n', (np.int64, n_frac)], ['vel', LagrangianVelocity], ['sigma', LagrangianVelocity]]
        if (calc_energy):
            keys = keys + [['epot',(np.float64, n_frac)], ['vr',(np.float64, n_frac)]]
            if (external_mode!='none'):
                keys = keys + [['epot_ext',(np.float64, n_frac)]]
        DictNpArrayMix.__init__(self, keys, _dat, _offset, _append, **kwargs)
        self.initargs['mass_fraction'] = m_frac
        self.initargs['calc_energy'] = calc_energy
        self.initargs['external_mode'] = external_mode
        self.initargs['calc_multi_rc'] = calc_multi_rc

        #keys = [['r', (np.float64, n_frac)],['m', (np.float64, n_frac)],['n', (np.float64, n_frac)]] # radius, mass, number
        #DictNpArrayMix.__init__(self, keys, _dat, _offset, _append, **kwargs)
        #self.vel  = LagrangianVelocity(_dat, _offset+self.ncols, False, **kwargs)
        #self.ncols += self.vel.ncols
        #self.keys.append(['vel',LagrangianVelocity])
        #self.sigma= LagrangianVelocity(_dat, _offset+self.ncols, False, **kwargs)
        #self.ncols += self.sigma.ncols
        #self.keys.append(['sigma',LagrangianVelocity])
        #self.initargs['mass_fraction'] = m_frac

    def calcTrh(self, G, gamma=0.02, mode='sphere'):
        """ Calculate Spitzer one-component half-mass relaxation time
            Trh = 0.138 N^0.5 Rh^1.5 /( G^0.5 m^0.5 ln(gamma N))
            Then add the result as a new member 'trh'.

            Notice mass fraction must contain 0.5 to use this function
            m use the average mass within half-mass radius

        Parameters
        ----------
        G: float
           Gravitational constant
        gamma: float (0.02 # Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037)
           The coefficient for Coulomb logarithm
        mode: string (sphere)
            sphere: calculate averaged properties from center to Lagrangian radii
            shell: calculate properties between two neighbor Lagrangian radii

        """
        rhindex=np.where(self.initargs['mass_fraction']==0.5)[0][0]
        if (mode=='shell'):
            n=self.n[:,0:(rhindex+1)].sum(axis=1)
            m=(self.m[:,0:(rhindex+1)]*self.n[:,0:(rhindex+1)]).sum(axis=1)/n
            trh=calcTrh(2*n, self.r[:,rhindex], m, G, gamma)
            self.addNewMember('trh',trh)
        else:
            trh=calcTrh(2*self.n[:,rhindex], self.r[:,rhindex], self.m[:,rhindex], G, gamma)
            self.addNewMember('trh',trh)

    def calcTcr(self, G, mode='sphere'):
        """ Calculate half-mass crossing time
            Tcr = Rh^1.5/sqrt(G M)
            Then add the result as a new member 'tcr'.

            Notice mass fraction must contain 0.5 to use this function

        Parameters
        ----------
        G: float
           Gravitational constant
        mode: string (sphere)
            sphere: calculate averaged properties from center to Lagrangian radii
            shell: calculate properties between two neighbor Lagrangian radii

        """
        rhindex=np.where(self.initargs['mass_fraction']==0.5)[0][0]
        if (mode=='shell'):
            M=(self.m[:,0:(rhindex+1)]*self.n[:,0:(rhindex+1)]).sum(axis=1)*2
            tcr=calcTcr(M, self.r[:,rhindex], G)
            self.addNewMember('tcr',tcr)
        else:
            tcr=calcTcr(self.n[:,rhindex]*self.m[:,rhindex]*2, self.r[:,rhindex], G)
            self.addNewMember('tcr',tcr)

    def calcOneSnapshot(self, _particle, _rc, _mode='sphere', **kwargs):
        """ Calculate one snapshot lagrangian parameters

        Parameters
        ----------
        _particle: inherited SimpleParticle
            particles data set sorted by distance to the coordinate center
        _rc: float
            core radius
        _mode: string (sphere)
            sphere: calculate averaged properties from center to Lagrangian radii
            shell: calculate properties between two neighbor Lagrangian radii
        ----------
        keyword arguments:
            read_rlagr: 1D np.array([])
                 When this option is used, Lagragian radii are not calculated,
                 but are given by the argument (an array that contain Lagragian radii).
                 Other properties (e.g. average mass, velocity, disperion) are calculated based on these.
        ----------
        return
            rlagr: 1D np.array
                Lagrangian radii of current data
        """
        shell_mode = True if (_mode == 'shell') else False

        mass_fraction = self.initargs['mass_fraction']
        calc_energy = self.initargs['calc_energy']
        external_mode = self.initargs['external_mode']
        calc_multi_rc = self.initargs['calc_multi_rc']
        n_frac = mass_fraction.size+1

        rc = _rc
        if calc_multi_rc:
            if (_particle.size>6):
                core = Core()
                kdtree = sp.cKDTree(_particle.pos)
                cm_pos, cm_vel = core.calcDensityAndCenter(_particle, kdtree)
                _particle.correctCenter(cm_pos, cm_vel)
                _particle.calcR2()
                rc = core.calcCoreRadius(_particle)
            else:
                rc = 0

        def find_mass_index(_mass_cum,mass_fraction):
            mass_bins=np.append(0,mass_fraction*_mass_cum[-1])
            index,count=np.histogram(_mass_cum,bins=mass_bins)
            index=index.cumsum()
            index[index>=_mass_cum.size]=_mass_cum.size-1
            return index

        if (_particle.size<=1):
            size = self.size
            empty = Lagrangian(np.zeros([1,n_frac*self.ncols]),**self.initargs)
            self.append(empty)
            if (size+1!=self.size):
                raise ValueError('Size should increase one, but increase', self.size-size)
            return np.zeros(mass_fraction.size)
        else:
            self.size += 1
            self.vel.size += 1
            self.sigma.size += 1
            r = np.sqrt(_particle.r2)
            rindex=np.array([])
            mcum=_particle.mass.cumsum()
            if ('read_rlagr' in kwargs.keys()):
                rbin = np.append(0,kwargs['read_rlagr'])
                rindex, e = np.histogram(r, bins=rbin)
                rindex=rindex.cumsum()
                rindex[rindex>=r.size]=r.size-1
            else:
                rindex= find_mass_index(mcum, mass_fraction)
            rlagr = r[rindex]
            if(len(self.r.shape)!=2):
                raise ValueError('r shape is wrong',self.r.shape)
            self.r = np.append(self.r, [np.append(rlagr,rc)], axis=0)
            nlagr = rindex+1
            if (shell_mode): nlagr[1:] -= nlagr[:-1]
            nc = (_particle.r2<(rc*rc)).sum()
            self.n = np.append(self.n, [np.append(nlagr,nc)], axis=0)
            mlagr = mcum[rindex]
            if (shell_mode): mlagr[1:] -= mlagr[:-1]
            sel = nlagr>0
            mlagr[sel] /= nlagr[sel]
            if (nc>0): mc = mcum[nc-1]/nc
            else:      mc = 0.0
            self.m = np.append(self.m, [np.append(mlagr,mc)], axis=0)

            m   =_particle.mass
            vel = _particle.vel
            pos = _particle.pos
            rx = pos[:,0]
            ry = pos[:,1]
            rz = pos[:,2]
            vx = vel[:,0]
            vy = vel[:,1]
            vz = vel[:,2]

            # x-y plane radial velocity * rxy
            rvxy = rx*vx + ry*vy
            # radial velocity value
            vr   = (rvxy + rz*vz)/r
            # tangential velocity vector
            vt = [None]*3
            vt[0] = vx - vr*rx/r
            vt[1] = vy - vr*ry/r
            vt[2] = vz - vr*rz/r
            # x-y plane radial position square
            rxy2 = rx*rx + ry*ry
            # rotational velocity
            vrotx = vx - rvxy*rx/rxy2
            vroty = vy - rvxy*ry/rxy2
            vrot = np.sqrt(vrotx*vrotx+vroty*vroty)
            # rotational direction sign
            vrotd = vrotx*ry - vroty*rx
            ng_sig = (vrotd<0.0)
            vrot[ng_sig] = -vrot[ng_sig]

            n_offset = np.append(np.zeros(1),nlagr).astype(int)
            vlst = [vx, vy, vz, vr, vt[0], vt[1], vt[2], vrot]
            vave = [None]*len(vlst)
            for k in range(len(vlst)):
                vlagr = []
                if (shell_mode):
                    vlagr = [np.average(m[n_offset[i]:n_offset[i+1]]*vlst[k][n_offset[i]:n_offset[i+1]])/mlagr[i] if (n_offset[i]<n_offset[i+1]) else 0.0 for i in range(mlagr.size)]
                else:
                    vlagr = [np.average(m[:n_offset[i+1]]*vlst[k][:n_offset[i+1]])/mlagr[i] if (n_offset[i+1]>0) else 0.0 for i in range(mlagr.size)]
                # core radius
                vlagr.append(np.average(m[0:nc]*vlst[k][0:nc])/mc if (nc>0) else 0.0)
                vave[k] = np.array(vlagr)

            self.vel.x = np.append(self.vel.x, [vave[0]], axis=0)
            self.vel.y = np.append(self.vel.y, [vave[1]], axis=0)
            self.vel.z = np.append(self.vel.z, [vave[2]], axis=0)
            self.vel.abs= np.append(self.vel.abs, [np.sqrt(vave[0]*vave[0]+vave[1]*vave[1]+vave[2]*vave[2])], axis=0)
            self.vel.rad = np.append(self.vel.rad, [vave[3]], axis=0)
            self.vel.tan = np.append(self.vel.tan, [np.sqrt(vave[4]*vave[4]+vave[5]*vave[5]+vave[6]*vave[6])], axis=0)
            self.vel.rot = np.append(self.vel.rot, [vave[7]], axis=0)

            sigma = [None]*len(vlst)
            for k in range(len(vlst)):
                slagr = []
                if (shell_mode):
                    slagr = [np.average(m[n_offset[i]:n_offset[i+1]] * (vlst[k][n_offset[i]:n_offset[i+1]] - vave[k][i])**2) / mlagr[i] if (n_offset[i]<n_offset[i+1]) else 0.0 for i in range(mass_fraction.size)]
                else:
                    slagr = [np.average(m[:n_offset[i+1]] * (vlst[k][:n_offset[i+1]] - vave[k][i])**2) / mlagr[i] if (n_offset[i+1]>0) else 0.0 for i in range(mass_fraction.size)]
                # core radius
                slagr.append(np.average(m[0:nc] * (vlst[k][0:nc] - vave[k][-1])**2) / mc if (nc>0) else 0.0)
                sigma[k] = np.array(slagr)

            sigma_all = sigma[0]+sigma[1]+sigma[2]
            self.sigma.x = np.append(self.sigma.x, [np.sqrt(sigma[0])], axis=0)
            self.sigma.y = np.append(self.sigma.y, [np.sqrt(sigma[1])], axis=0)
            self.sigma.z = np.append(self.sigma.z, [np.sqrt(sigma[2])], axis=0)
            self.sigma.abs= np.append(self.sigma.abs, [np.sqrt(sigma_all)], axis=0)
            self.sigma.rad = np.append(self.sigma.rad, [np.sqrt(sigma[3])], axis=0)
            self.sigma.tan = np.append(self.sigma.tan, [np.sqrt(sigma[4]+sigma[5]+sigma[6])], axis=0)
            self.sigma.rot = np.append(self.sigma.rot, [np.sqrt(sigma[7])], axis=0)

            if (calc_energy):
                epot = []
                if (shell_mode):
                    epot = [np.sum(_particle.pot[n_offset[i]:n_offset[i+1]] * _particle.mass[n_offset[i]:n_offset[i+1]]) if (n_offset[i]<n_offset[i+1]) else 0.0 for i in range(mass_fraction.size)]

                else:
                    epot = [np.sum(_particle.pot[:n_offset[i+1]] * _particle.mass[:n_offset[i+1]]) if (n_offset[i+1]>0) else 0.0 for i in range(mass_fraction.size)]

                epot.append(np.sum(_particle.pot[:nc]*_particle.mass[:nc]  if (nc>0) else 0.0))
                ekin = sigma_all*np.append(mlagr,mc)*np.append(nlagr,nc)
                epot = np.array(epot)
                self.epot = np.append(self.epot, [epot], axis=0)
                self.vr = np.append(self.vr, [-2.0*ekin/epot], axis=0)
                if (external_mode != 'none'):
                    epot_ext = []
                    if (shell_mode):
                        epot_ext = [np.sum(_particle.pot_ext[n_offset[i]:n_offset[i+1]] * _particle.mass[n_offset[i]:n_offset[i+1]]) if (n_offset[i]<n_offset[i+1]) else 0.0 for i in range(mass_fraction.size)]
                    else:
                        epot_ext = [np.sum(_particle.pot_ext[:n_offset[i+1]] * _particle.mass[:n_offset[i+1]]) if (n_offset[i+1]>0) else 0.0 for i in range(mass_fraction.size)]
                    epot_ext.append(np.sum(_particle.pot_ext[:nc]*_particle.mass[:nc]  if (nc>0) else 0.0))
                    self.epot_ext = np.append(self.epot_ext, [np.array(epot_ext)], axis=0)
            return rlagr

class LagrangianMultiple(DictNpArrayMix):
    """ Lagrangian for single, binaries and all
    Keys: (class members)
        single (Lagrangian): Lagrangian data for single particles
        binary (Lagrangian): Lagrangian data for binaries
        all (Lagrangian): Lagrangian data for all data (binary is treated as c.m., count once)
        Additional members that depend on the keyword argument 'add_star_type':
            [add_star_type] (Lagrangian): Lagrangian data for specific type of stars, see the following content for details

    Keyword argument 'add_star_type' (list of string):
        In case when interrupt_mode='(mo)bse', a list of star types can be given,
        so that the Lagrangian properties for these specific types will be calculated.
        There are four styles of star types:
           (1) a single SSE type name (see help(petar.SSEType))
               For example, if 'BH' is given, one additional class member (type is class Lagrangian), BH, is added.
               The Lagrangian radii are determined by only counting BHs,
               then the properties (average mass, velocity...) of BHs are calculated using these Lagrangian radii.
           (2) a combination of different SSE types connected by '_'.
               This will include mutliple SSE types in the calculation.
               For example, if 'BH_NS_WD' is given, the additional class member, BH_NS_WD,
               count BHs, NSs and WDs together to calculate the Lagrangian properties.
           (3) a single SSE type name with the prefix 'no'
               This will exclude the given SSE type name in the calculation.
               For exmaple, if 'noBH' is given, the additional class member, noBH,
               count all stars except BHs to calculate the Lagrangian properties.
           (4) two types (can be any case of the style 1-3) are given by '[type 1]__in__[type 2]'
               Then, the properties of type 1 stars within the spheres or shells of
               Lagragian radii of type 2 stars are calculated. The Lagrangian radii of type 1 are not calculated.
               For example, if 'BH__in__all' is given, Lagrangian properties of BHs are calculated
               within the shell or sphere of Lagrangian radii of all stars (instead of Lagrangian radii of BHs).
               Notice that the two type names (except 'all') should also be added separately in the list.
               In the case of 'BH__in__all', add_star_type must contain 'BH', i.e. add_star_type=['BH','BH__in__all', ...].
               In another example, 'BH__in__MS', add_star_type=['BH','MS','BH__in__MS',...].

    For binaries:
        One binary is counted once in class member 'binary', its c.m. data is used.

        In case of [add_star_type], one binary is counted once using its c.m. position and velocity.
        If two components are both this type of star, the total mass is used.
        If one component is this type of star, its mass and c.m. position and velocity are used.

    """
    def __init__ (self, _dat=None, _offset=int(0), _append=False, **kwargs):
        """ DictNpArrayMix type initialzation, see help(DictNpArrayMix.__init__)

        Parameters
        ----------
        keyword arguments:
            mass_fraction: 1D numpy.ndarray (np.array([0.1, 0.3, 0.5, 0.7, 0.9]))
                An array to indicate the mass fractions to calculate lagrangian radii.
            calc_energy: bool (False)
                If true, add kinetic, potential energies and virial ratio calculation
            external_mode: string (none)
                PeTar external mode (set in configure): galpy, none
                If it is not none, epot_ext will be added if calc_energy=True
            calc_multi_rc: bool (False)
                If true, calculate core radius using KDTree and correct the central position
            add_star_type: list of string ([])
                An array containing the star type names to calculate additional Lagrangian properties for specific star types.
        """
        m_frac=np.array([0.1,0.3,0.5,0.7,0.9])
        if ('mass_fraction' in kwargs.keys()): m_frac=kwargs['mass_fraction'].copy()
        calc_energy=False
        if ('calc_energy' in kwargs.keys()): calc_energy=kwargs['calc_energy']
        external_mode='none'
        if ('external_mode' in kwargs.keys()): external_mode=kwargs['external_mode']
        calc_multi_rc=False
        if ('calc_multi_rc' in kwargs.keys()): calc_multi_rc=kwargs['calc_multi_rc']

        keys=[['time',np.float64], ['single',Lagrangian], ['binary', Lagrangian], ['all', Lagrangian]]
        add_star_type=[]
        if ('add_star_type' in kwargs.keys()):
            add_star_type=kwargs['add_star_type'].copy()
            for name in add_star_type:
                keys += [[name, Lagrangian]]

        DictNpArrayMix.__init__(self, keys, _dat, _offset, _append, **kwargs)
        self.initargs['mass_fraction'] = m_frac
        self.initargs['calc_energy'] = calc_energy
        self.initargs['external_mode'] = external_mode
        self.initargs['calc_multi_rc'] = calc_multi_rc
        self.initargs['add_star_type'] = add_star_type

        #DictNpArrayMix.__init__(self, [['time',np.float64]], _dat, _offset, _append, **kwargs)
        #self.single = Lagrangian(_dat, _offset+self.ncols, False, **kwargs)
        #self.ncols += self.single.ncols
        #self.keys.append(['single',Lagrangian])
        #self.binary = Lagrangian(_dat, _offset+self.ncols, False, **kwargs)
        #self.ncols += self.binary.ncols
        #self.keys.append(['binary',Lagrangian])
        #self.all    = Lagrangian(_dat, _offset+self.ncols, False, **kwargs)
        #self.ncols += self.all.ncols
        #self.keys.append(['all',Lagrangian])
        #self.size = self.all.size

    def calcOneSnapshot(self, time, single, binary, rc, mode):
        """ Calculate Lagrangian radii and related properties for one snapshot

        Parameters
        ----------
        time: float
            current evolved time of the system
        single: inherited SimpleParticle
            single partilces (center is corrected and r2 exist)
        binary: Binary
            binaries (center is corrected and r2 exist)
        rc: float
            Core radius
        mode: string
            sphere: calculate averaged properties from center to Lagrangian radii
            shell: calculate properties between two neighbor Lagrangian radii
        """
        self.time = np.append(self.time, time)
        single_sim = SimpleParticle(single)
        single_sim.calcR2()
        single_sim.addNewMember('pot',single.pot)
        if (self.initargs['external_mode']!='none'):
            single_sim.addNewMember('pot_ext',single.pot_ext)
        binary_sim = SimpleParticle(binary)
        binary_sim.calcR2()
        binary.calcPot()
        binary_sim.addNewMember('pot',binary.pot)
        if (self.initargs['external_mode']!='none'):
            binary.calcPotExt()
            binary_sim.addNewMember('pot_ext',binary.pot_ext)
        if len(self.initargs['add_star_type'])>0:
            single_sim.addNewMember('type1',single.star.type)
            single_sim.addNewMember('m1',single.mass)
            single_sim.addNewMember('type2',single.star.type)
            single_sim.addNewMember('m2',np.zeros(single.size))
            binary_sim.addNewMember('type1',binary.p1.star.type)
            binary_sim.addNewMember('m1',binary.p1.mass)
            binary_sim.addNewMember('type2',binary.p2.star.type)
            binary_sim.addNewMember('m2',binary.p2.mass)
        all_sim = join(single_sim, binary_sim)
        n_single = single.size
        n_binary = binary.size

        idx = all_sim.r2.argsort()
        idx_single = idx[idx<n_single]
        idx_binary = idx[idx>=n_single] - n_single
        single_sort = single_sim[idx_single]
        binary_sort = binary_sim[idx_binary]
        all_sort = all_sim[idx]

        self.single.calcOneSnapshot(single_sort, rc, mode)
        self.binary.calcOneSnapshot(binary_sort, rc, mode)
        rlagr=dict()
        rlagr['all'] = self.all.calcOneSnapshot(all_sort, rc, mode)
        if (self.binary.size != self.single.size):
            raise ValueError('Size inconsistence: single.size:', self.single.size, ' binary.size:', self.binary.size)

        name_list=[]
        cross_list=[]
        obj_sort=dict()
        for name_org in self.initargs['add_star_type']:
            name_cross= name_org.split('__in__')
            if (len(name_cross)>1):
                cross_list.append(name_cross)
                obj_sort[name_cross[0]]=np.array([])
            else:
                name_list.append(name_cross[0])

        for name in name_list:
            sel = np.zeros(all_sort.size).astype(bool)
            all_sort.mass = np.zeros(all_sort.size) # clear up mass, final mass is the sum of matched component masses
            for subname in name.split('_'):
                sel1 = np.array([]).astype(bool)
                sel2 = np.array([]).astype(bool)
                if (subname[:2]=='no'):
                    type_index = BSE_STAR_TYPE_INDEX[subname[2:]]
                    sel1 = (all_sort.type1!=type_index)
                    sel2 = (all_sort.type2!=type_index)
                else:
                    type_index = BSE_STAR_TYPE_INDEX[subname]
                    sel1 = (all_sort.type1==type_index)
                    sel2 = (all_sort.type2==type_index)
                sel = sel | (sel1 | sel2)
                all_sort.mass[sel1] += all_sort.m1[sel1]
                all_sort.mass[sel2] += all_sort.m2[sel2]
            all_sel_sort = all_sort[sel]
            rlagr[name]=self.__dict__[name].calcOneSnapshot(all_sel_sort, rc, mode)
            # save particle list for cross check
            if (name in obj_sort.keys()): obj_sort[name] = all_sel_sort
            #self.__dict__[name+'_cross'].calcOneSnapshot(all_sel_sort, rc, mode, read_rlagr=rlagr)

        # cross check
        for name_pair in cross_list:
            self.__dict__[name_pair[0]+'__in__'+name_pair[1]].calcOneSnapshot(obj_sort[name_pair[0]], rc, mode, read_rlagr=rlagr[name_pair[1]])

        self.size += 1

max_stars_repo_name: michelamapelli/PeTar

max_issues_repo_name: michelamapelli/PeTar

max_forks_repo_name: michelamapelli/PeTar

reason: import numpy,from scipy,import astropy,from astropy


text:
import numpy as np
import matplotlib.pyplot as plt

### numerical model specifications and physical constants
def getPhysConstants():
    physConst = {}
    physConst['RE'] = 6371.0             # Earth's radius in km
    physConst['sigSB'] = 5.6703*(10**-8) # Stefan-Boltzmann constant
    physConst['KelToCel'] = 273.15       # 0 C in Kelvin
    return physConst

def getInputParamsDefault():
    params = {}
    physConst = getPhysConstants()

    ## environment (radiation from Sun and Earth)
    params['Fsun'] = 1372   # solar flux, W/m2
    params['FIR'] = 240     # Earth's IR flux, W/m2
    params['rhoE'] = 0.3    # Earth's reflectivity

    ## orbit
    params['h'] = 550.0     # in km, orbit's altitude
    params['PorbMin'] = 90  # orbital period in minutes
    params['etaP'] = 0.33   # fraction of orbit in eclipse
    params['falb'] = 0.62   # correction for albedo variation
                            # for polar orbit facing Sun: 0.06
    # view factor
    RE = physConst['RE']
    params['fE'] = (RE/(RE+params['h']))**2

    ## satellite
    # area
    params['Atot'] = 0.10   # in m2, as 2U CubeSat
    params['etaS'] = 0.25   # effective area ratio, sphere, point source
    params['etaE'] = 0.50   # effective area ratio, sphere, hemisphere
    # use instead the version with a correction for incomplete hemisphere
    # it's basically viewing factor for two spheres multiplied by fE
    params['etaE'] = 0.50 * (1-np.sqrt(1-params['fE'])) / params['fE']

    # surface emissivity
    params['alphaS'] = 0.86 # black anodized Al
    params['epsT'] = 0.86   # black anodized Al
    params['alphaIR'] = params['epsT']  # usual ansatz
    # thermal inertia
    params['mass'] = 2.0    # kg (n.b. implies density)
    params['C'] = 921.0     # J/kg/K, aluminum
    # battery charging/power dissipation
    params['etaCell'] = 0.2 # fraction of energy for charging

    return params

def getInputParamsHot(params=""):

    if (params==""):
        params = getInputParamsDefault()

    ## modify for the hot case
    # environment (radiation from Sun and Earth)
    params['Fsun'] = 1422   # solar flux, W/m2
    params['FIR'] = 260     # Earth's IR flux, W/m2
    params['rhoE'] = 0.35   # Earth's reflectivity
    ## orbit
    params['etaP'] = 0.0    # no eclipse in beta=90 (terminator) orbit
    params['falb'] = 0.06   # correction for albedo in hot orbit

    return params

def getInputParamsCold(params=""):

    if (params==""):
        params = getInputParamsDefault()

    ## modify for the cold case
    # environment (radiation from Sun and Earth)
    params['Fsun'] = 1322   # solar flux, W/m2
    params['FIR'] = 220     # Earth's IR flux, W/m2
    params['rhoE'] = 0.25   # Earth's reflectivity

    return params

### heat flux and equilibrium temperature evaluations

def getQsun(p):
    return p['etaS']*p['Atot']*p['alphaS']*p['Fsun']

def getQref(p):
    return p['fE']*p['etaE']*p['Atot']*p['falb']*p['rhoE']*p['alphaS']*p['Fsun']

def getQIR(p):
    return p['fE']*p['etaE']*p['Atot']* p['alphaIR']*p['FIR']

def getQdissip(p):
    return p['etaCell']*(1-p['etaP'])*(getQsun(p)+getQref(p))

def getQin(p):
    QinSun = (1-p['etaCell']*p['etaP'])*(getQsun(p)+getQref(p)) + getQIR(p)
    QinEclipse = getQIR(p) + getQdissip(p)
    return QinSun, QinEclipse

def getAllHeatQ(p):
    return getQsun(p), getQref(p), getQIR(p), getQdissip(p)

## equilibrium temperatures
def getTeq(Qin,p):
    physConst = getPhysConstants()
    return (Qin/(p['Atot']*p['epsT']*physConst['sigSB']))**0.25

# return the view factor divided by fE=(RE/(RE+altitude))^2
# here fE = 1/h^2
def getF12h2(h,beta):
    rad2deg = 180/np.pi
    brad = beta/rad2deg
    if (brad<=np.arccos(1/h)):
        return np.cos(brad)
    x = np.sqrt(h**2-1)
    y = -x / np.tan(brad)
    if (y>1):
        return 0
    z = np.sqrt(1-y**2)
    t1 = (np.cos(brad)*np.arccos(y)-x*np.sin(brad)*z)
    t2 = np.arctan(np.sin(brad)*z/x)
    # print(x,y,z,t1,t2)
    return (t1+t2*h**2)/np.pi

def getEffectiveAreas2U(h, beta):
    if ((beta<0)or(beta>180)):
        return 0
    if (beta>90):
        beta = 180-beta
    # small side towards Earth
    s1 = getF12h2(h, beta)
    # the "other" small side
    s2 = getF12h2(h, 180-beta)
    # long 2U side "towards" Earth
    s3 = 2 * getF12h2(h, beta+90)
    # long 2U side "away from" Earth
    s4 = 2 * getF12h2(h, 180 - (beta+90))
    # 2 long 2U sides "perpendicular to" Earth
    s5 = 2 * 2 * getF12h2(h, 90)
    # print(h,beta,s1, s2, s3, s4, s5)
    return (s1+s2+s3+s4+s5)/10

def getEffectiveAreas1U(h, beta):
    if ((beta<0)or(beta>180)):
        return 0
    if (beta>90):
        beta = 180-beta
    # small side towards Earth
    s1 = getF12h2(h, beta)
    # the "other" small side
    s2 = getF12h2(h, 180-beta)
    # small side "towards" Earth
    s3 = getF12h2(h, beta+90)
    # small side "away from" Earth
    s4 = getF12h2(h, 180 - (beta+90))
    # 2 small sides "perpendicular to" Earth
    s5 = 2 * getF12h2(h, 90)
    # print(h,beta,s1, s2, s3, s4, s5)
    return (s1+s2+s3+s4+s5)/6

## for compatibility with old code

def doOneCase(modelTitleText,alpha,epsilon, params="", verbose=False):
    physConst = getPhysConstants()
    titlePrint(modelTitleText)
    if (params == ""):
        print('using default parameters from getInputParamsDefault()')
        params = getInputParamsDefault()
    params['alphaS'] = alpha
    params['epsT'] = epsilon
    params['alphaIR'] = params['epsT']

    # heat fluxes
    Q_sun_flux, Q_albedo_flux, Q_IR_flux, Q_int = getAllHeatQ(params)
    if (verbose):
        print('Q (sun, alb, IR, diss):', Q_sun_flux, Q_albedo_flux, Q_IR_flux, Q_int)

    # equilibrium solutions
    Q_hot, Q_cold = getQin(params)
    print('Qsun=', Q_hot, ' Qeclipse=', Q_cold)
    Temp_hot = getTeq(Q_hot,params)
    Temp_cold = getTeq(Q_cold,params)
    if (verbose):
        printEqTemp('Hot', Temp_hot, physConst['KelToCel'])
        printEqTemp('Cold', Temp_cold, physConst['KelToCel'])

    # solve for time variation
    t0H = getTimeConstant(params,Temp_hot)
    t0C = getTimeConstant(params,Temp_cold)
    xH = (1-params['etaP'])*params['PorbMin']*60/t0H
    xC = params['etaP']*params['PorbMin']*60/t0C
    Tmin, Tmax, tau0C, tauFC, tau0H, tauFH = solveBistableEquilibrium(xC, xH, Temp_cold, Temp_hot, 1000)

    # generate temperature arrays
    timeA, TempsA = getFullTempArray(params, Tmin, Tmax, Temp_cold, Temp_hot, t0C, t0H)
    if (verbose):
        print('doOneCase: temp. range =', np.min(TempsA), np.max(TempsA))
        print('     in Celsius: range =', np.min(TempsA-273.15), np.max(TempsA-273.15))

    return timeA, TempsA

def getTimeConstant(params,Teq):
    physConst = getPhysConstants()
    return (params['mass']*params['C']/physConst['sigSB']/params['Atot']/params['epsT']/Teq**3)

def getFullTempArray(params, Tmin, Tmax, Temp_cold, Temp_hot, t0C, t0H):

    # generate temperature arrays
    if (Tmin >= Temp_hot):
        print('adjusting Tmin=', Tmin, 'to: Thot=', Temp_hot)
        Tmin = 0.998*Temp_hot
    if (Tmin <= Temp_cold):
        print('adjusting Tmin=', Tmin, 'to: Tcold=', Temp_cold)
        Tmin = 1.001*Temp_cold
    if (Tmax >= Temp_hot):
        print('adjusting Tmax=', Tmax, 'to: Thot=', Temp_hot)
        Tmax = 0.999*Temp_hot
    TempsC = np.linspace(Tmax, Tmin, 100)
    timeC = getTempSolution(TempsC, t0C, Temp_cold)
    TempsH = np.linspace(Tmin, Tmax, 200)
    timeH = getTempSolution(TempsH, t0H, Temp_hot)
    timeH[-1] = (1-params['etaP'])*params['PorbMin']*60

    # concatenate both solutions
    timeOffset = params['etaP']*params['PorbMin']*60
    timeHot = (1-params['etaP'])*params['PorbMin']*60
    timeHshifted = timeH/np.max(timeH)*timeHot + timeOffset
    timeA = np.concatenate((timeC/np.max(timeC)*timeOffset, timeHshifted), axis=None)
    TempsA = np.concatenate((TempsC, TempsH), axis=None)

    return timeA, TempsA

def getNumSolution(params, Tstart, Tc=0.0, Pc=0.0):

    orbPeriodMin = params['PorbMin']
    etaP = params['etaP']
    t_final_min = orbPeriodMin*etaP   # eclipse duration in minutes
    t_step = 1.0                  # integration step in seconds

    # Tstart = TempsBlack[0]
    Qs, Qe = getQin(params)
    epsT = params['epsT']
    Atot = params['Atot']
    C = params['C']
    mass = params['mass']

    physConst = getPhysConstants()
    sigmaSB = physConst['sigSB']

    if (etaP>0):
        facSun = (1-etaP)/etaP
    else:
        facSun = 1.0

    # need to enforce cyclic boundary condition with iterations
    # the duration of "sunshine", in units of eclipse time (=t_final_min)
    for k in range(0,25):
        time1, Temps1 = nonSteadyStateTemp(Tstart,t_final_min,t_step,Qe,C,mass,sigmaSB,epsT,Atot,Tc,Pc)
        tmax = facSun*t_final_min
        time2, Temps2 = nonSteadyStateTemp(Temps1[-1],tmax,t_step,Qs,C,mass,sigmaSB,epsT,Atot,Tc,Pc)
        Tstart = 0.5*(Tstart+Temps2[-1])

    time = np.concatenate((time1, (time2+time1[-1])), axis=None)
    Temperature = np.concatenate((Temps1, Temps2), axis=None)

    return time, Temperature

# given x and a tau0 array, find tau for each tau0
def getTauSolution(x, tau0, Nsteps=100):
    tauFinal = 0*tau0
    for i in range(0,np.size(tau0)):
        thisTau0 = tau0[i]
        tauGrid = np.linspace(thisTau0, 1, Nsteps)
        notdone = True
        for j in range(1,np.size(tauGrid)):
            if (notdone):
                thisTau = tauGrid[j]
                thisX = getXforTempSolution(thisTau, thisTau0)
                if (thisX >= x):
                    tauFinal[i] = thisTau
                    notdone = False
    return tauFinal

# given tau and tau0, return x=t/t0
def getXforTempSolution(tau, tau0):
    x1 = 0.5*(np.arctan(tau)-np.arctan(tau0))
    x2 = 0.25*(np.log((tau+1)/(tau0+1)) - np.log((tau-1)/(tau0-1)))
    return x1+x2

# given an array of temperatures, Teq and t0, return model time grid
def getTempSolution(T, t0, Teq):
    tau = T/Teq
    # print('getTempSolution: Teq, t0, tau0, tauFinal=', Teq, t0, tau[0], tau[-1])
    t1 = 0.5*(np.arctan(tau)-np.arctan(tau[0]))
    t2 = 0.25*(np.log((tau+1)/(tau[0]+1)) - np.log((tau-1)/(tau[0]-1)))
    return t0*(t1+t2)

def solveBistableEquilibrium(xC, xH, TeqCold, TeqHot, Nsteps=100):
    C1 = TeqCold / TeqHot
    tau0Carr = np.linspace(1.01,3.0, Nsteps)
    tauFinalCarr = getTauSolution(xC, tau0Carr)
    tau0Harr = C1*tauFinalCarr
    tauFinalHarr = getTauSolution(xH, tau0Harr)
    tau0C2 = tauFinalHarr/C1
    ## now solve for tauFinalCarr(tau0Carr) = tauFinalCarr(tau0C2)
    ## to get tau0C and tauFinalC, and then tau0H, tauFinalH, and
    ## finally TminEq and TmaxEq
    # by construction, tau0Carr and tau0C2 are increasing
    # and at first point tau0C2[0] > tau0Carr[0]
    notdone = True
    for i in range(1,np.size(tau0C2)):
        if (notdone):
            if (tau0Carr[i] >= tau0C2[i]):
                tau0C = tau0Carr[i]
                tauFinalC = tauFinalCarr[i]
                notdone = False
    tau0H = C1 * tauFinalC
    tauFinalH = C1 * tau0C
    TminEq = tauFinalC * TeqCold
    TmaxEq = tau0C * TeqCold
    return TminEq, TmaxEq, tau0C, tauFinalC, tau0H, tauFinalH

############# old code from Haley Stewart (University of Washington)   #################################
def nonSteadyStateTemp(Tinitial,t_final_min,t_step,Qheating,c,mass,sigma,emissivity,Area_sphere,Tc,Pc):
    T = Tinitial # set the initial temp to the hot eq temp
    Temps = []
    Temps.append(T) # include initial temp in list
    time = []
    time.append(0) # python counts from zero
    # need smaller time step
    t_final = t_final_min*60 # convert minutes to seconds
    for i in range(1, int((t_final/t_step) + 1), 1): # sampling every t_step seconds
        T = nextTemp(T, t_step, Qheating, c, mass, sigma, emissivity, Area_sphere,Tc,Pc)
        Temps.append(T)
        time.append(i*t_step)
    return np.array(time), np.array(Temps)

def nextTemp(T, t_step, Qheating, c, mass, sigma, emissivity, Area_sphere,Tc,Pc):
    if ((T<Tc)):
        QheatingTotal = Qheating + Pc
    else:
        QheatingTotal = Qheating
    dTdt = (1/(c*mass))*(QheatingTotal-(sigma*(T**4)*emissivity*Area_sphere))
    return T + dTdt*t_step
#########################################################################################################

### PLOTTING
def TempsPlot2(time1, temp1, c1, time2, temp2, c2, time3, temp3, c3, outfile="", title=""):
    physConst = getPhysConstants()
    KelToCel = physConst['KelToCel']

    # Plot the Temperature wrt time
    battTmin =  0.0
    battTmax = 40.0
    # operating range for ClydeSpace batteries
    battTmin = -10.0
    battTmax =  40.0
    tempmax = []
    tempmin = []
    for int in range(0,len(time1)):
        tempmax.append(battTmax)
        tempmin.append(battTmin)
    timeMin = np.array(time1/60)

    fig, ax = plt.subplots()
    fig.subplots_adjust(wspace=0.22, left=0.16, right=0.97, bottom=0.15, top=0.97)

    ax.tick_params(axis='both', which='major', labelsize=15)
    plt.plot(np.array(time1)/60,np.array(temp1)-KelToCel,label=c1, lw=3)
    plt.plot(np.array(time2)/60,np.array(temp2)-KelToCel,label=c2, lw=3)
    if (c3 != ""):
        plt.plot(np.array(time3)/60,np.array(temp3)-KelToCel,label=c3, lw=3)
    plt.plot(timeMin, battTmin+0*timeMin, lw=1, c='black')
    plt.plot(timeMin, battTmax+0*timeMin, lw=1, c='black')
    # ax.fill_between(timeMin, battTmin+0*timeMin, battTmax+0*timeMin, alpha=0.1)
    tt = np.linspace(0,96,960)
    ax.fill_between(tt, battTmin+0*tt, battTmax+0*tt, alpha=0.1)

    plt.rc('legend',fontsize=12)
    plt.legend(loc=1)
    plt.title(title)
    plt.xlabel('Time (minute)', fontsize=18)
    plt.ylabel('Temperature (C)', fontsize=18)
    yLimB, yLimT = plt.ylim()
    plt.ylim(yLimB-5, 60.0)

    if (outfile == ""):
        name = 'figures/TempsVsOperatingTemp.png'
    else:
        name = 'figures/TempsPlot_' + outfile + '.png'
    plt.savefig(name)
    plt.close("all")
    return

def TempsPlotCompare(timeA, TempsA, timeN, TempsN, labelText, label2, outfile="", title=""):
    physConst = getPhysConstants()
    KelToCel = physConst['KelToCel']

    # Plot the Temperature wrt time
    battTmin =  0.0
    battTmax = 40.0
    # operating range for ClydeSpace batteries
    battTmin = -10.0
    battTmax =  40.0
    tempmax = []
    tempmin = []
    for int in range(0,len(timeA)):
        tempmax.append(battTmax)
        tempmin.append(battTmin)
    timeMin = np.array(timeA)/60

    fig, ax = plt.subplots()
    fig.subplots_adjust(wspace=0.22, left=0.18, right=0.97, bottom=0.17, top=0.97)

    ax.tick_params(axis='both', which='major', labelsize=20)
    plt.plot(np.array(timeA)/60,np.array(TempsA)-KelToCel,label=labelText, lw=3, c='b')
    plt.plot(np.array(timeN)/60,np.array(TempsN)-KelToCel,label=label2, lw=2, c='r', linestyle='dashed')
    plt.plot(np.array(timeN)/60, battTmin+0*timeMin, lw=1, c='black')
    plt.plot(np.array(timeN)/60, battTmax+0*timeMin, lw=1, c='black')
    # ax.fill_between(timeMin, tempmin, tempmax, alpha=0.1)
    ax.fill_between(np.array(timeN)/60, battTmin, battTmax, alpha=0.1)

    plt.rc('legend',fontsize=15)
    plt.legend(loc=1)
    plt.title(title)
    plt.xlabel('Time (minute)', fontsize=22)
    plt.ylabel('Temperature (C)', fontsize=22)
    yLimB, yLimT = plt.ylim()
    plt.ylim(yLimB-5, 60.0)

    if (outfile == ""):
        name = 'figures/TempsVsOperatingTemp.png'
    else:
        name = 'figures/TempsPlotCompare_' + outfile + '.png'

    plt.savefig(name)
    plt.close("all")
    return

def TempsPlot(time, temp1, coating1, temp2, coating2, temp3, coating3, outfile="", title=""):
    TempsPlot2(time, temp1, coating1, time, temp2, coating2, time, temp3, coating3, outfile, title)
    return

def titlePrint(coating):
    print()
    print(coating + ' Sphere, Temperature Extremes')
    print('-------------------------------------------------------')

def printEqTemp(case, Temp, KelToCel):
    print(case +' Eq. Temp: ' + str('%.2f'%(Temp)) + 'K = ' + str('%.2f'%(Temp-KelToCel)) + u'\u2103')

##########################################################################################
#### additional 2nd generation tools for supporting variable satellite orientation
#### along beta=0 or 90 deg. sun-synchronous orbits

### 1) for computing effective areas and other factors along an orbit

def getAreasForOrbitImagingMode(XYphiAngle, beta0=True, nU=2, Norbit=3600):
    ## assume that -Z side with the camera is always pointing towards Earth
    ## and two possible sun-synchronous orbits: beta=0 and beta=90
    # aux helper
    def phiEaux(phiS, alpha):
        return 180-phiS - (1-phiS/90)*alpha

    # angle along the orbit
    alpha = np.linspace(0, 360, Norbit)
    # orbital height parameter (h = (RE + altitude)/RE)
    h = 1.0865
    if beta0 == True:
        print('assuming beta=0 sun-synchronous orbit')
        ### define angles towards the Sun for a beta=0 orbit
        if False:
            # this is for Z- always pointing towards nadir: imaging mode
            thetaS = np.where(alpha<180, alpha, 360-alpha)
            # angles towards Earth when Z- always pointing towards nadir
            thetaE = 180 + 0*alpha
            phiE = XYphiAngle + 0*alpha
        else:
            # this is for solar sensor always pointing towards the Sun: sun seeker
            thetaS = 90 + 0*alpha
            # thetaE varies
            thetaE = np.where(alpha>90, np.abs(alpha-270), 90+alpha)
            phiE = np.where(alpha<180, phiEaux(50,alpha), phiEaux(50,360-alpha))
        # always pointing towards the Sun at fixed input angle
        phiS = XYphiAngle + 0*alpha  # sun sensor on Y+ face
        # maximum angle for eclipse (per Gilmore, eq.2.7 and 550 km orbit):
        aMax = 66.96 # corresponds to 37% of orbit in eclipse, 36 min for 96 min orbit
        # aMax = 56.52 # corresponds to 31% of orbit in eclipse, for beta=45 orbit
        eclipseFactor = np.where(np.abs(180-alpha)<aMax, 0, 1)
        albedoFactor = np.where(alpha<180, GilmoreAlbedo(alpha), GilmoreAlbedo(360-alpha))
    else:
        print('assuming beta=90 sun-synchronous orbit')
        # define angles towards the Sun for a beta=90 orbit
        thetaS = 90 + 0*alpha
        phiS = XYphiAngle + 0*alpha
        # angles towards Earth same for both cases
        thetaE = 180 + 0*alpha
        phiE = XYphiAngle + 0*alpha
        # these are trivial in this case
        eclipseFactor = 1 + 0*alpha
        albedoFactor = GilmoreAlbedo(90) + 0*alpha

    # compute area for each of the six sides
    areaS = getEffectiveAreasCubeSat(1000, phiS, thetaS, nU)
    areaE = getEffectiveAreasCubeSat(h, phiE, thetaE, nU)
    # and effective areas
    etaS = getEffectiveEtasCubeSat(areaS, nU)
    etaE = getEffectiveEtasCubeSat(areaE, nU)

    return alpha, eclipseFactor, albedoFactor, areaS, areaE, etaS, etaE

def GilmoreAlbedo(thetaDeg):
    # from Gilmore (2002):
    return np.where(thetaDeg<99.99, (np.cos(0.9*thetaDeg/180*np.pi))**1.5, 0)

## given a direction defined by spherical coordinate angles phideg and thetadeg,
## and a CubeSat aligned with the XYZ frame, return the six angles between the normal
## direction at each face and this input direction; in other words, get the projection
## area factor (to be multiplied by the geometric area and absorptivity)
def getNormalAngles(phideg, thetadeg):
    rad2deg = 180/np.pi
    Cphi = np.cos(phideg/rad2deg)
    Sphi = np.sin(phideg/rad2deg)
    Ctheta = np.cos(thetadeg/rad2deg)
    Stheta = np.sin(thetadeg/rad2deg)
    beta = {}
    beta['Xp'] = np.arccos(Stheta*Cphi)*rad2deg
    beta['Xm'] = 180-beta['Xp']
    beta['Yp'] = np.arccos(Stheta*Sphi)*rad2deg
    beta['Ym'] = 180-beta['Yp']
    beta['Zp'] = np.arccos(Ctheta)*rad2deg
    beta['Zm'] = 180-beta['Zp']
    return beta

## given the variation of phi and theta angles along an orbit,
## return the effective areas for each side (assumes a string
## of nU basic 1U units)
def getEffectiveAreasCubeSat(h, phideg, thetadeg, nU=2):
    angleArr = getNormalAngles(phideg, thetadeg)
    area = {}  # in units of 1U areas
    for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
        area[s] = getArea(h, angleArr, s)
    for s in ['Xp', 'Xm', 'Yp', 'Ym']:
        area[s] *= nU
    return area

def getEffectiveEtasCubeSat(area, nU=2):
    eta = area['Zp'] + area['Zm']
    totArea = 2
    for s in ['Xp', 'Xm', 'Yp', 'Ym']:
        eta += area[s]
        totArea += nU
    return eta/totArea

def getArea(h, angles, side):
    betaArr = angles[side]
    Nsteps = np.size(betaArr)
    A = np.zeros(Nsteps)
    for i in range(0,Nsteps):
        beta = betaArr[i]
        A[i] = getF12h2(h, beta)
    return A

def getEclipseFraction(betaOrbit, hOrbit=550):
    # Gilmore (2002), eq. 2.7
    physConst = getPhysConstants()
    RE = physConst['RE'] # in km, so hOrbit in km, too
    betaStar = np.arcsin(RE/(RE+hOrbit))/np.pi*180
    if betaOrbit > betaStar:
        return 0
    else:
        arg = np.sqrt(hOrbit**2+2*RE*hOrbit)/(RE+hOrbit)/np.cos(betaOrbit/180*np.pi)
        fEclipse = np.arccos(arg)/np.pi
    return fEclipse

### 2) for computing heating sources along an orbit

# note that area returned by getEffectiveAreasCubeSat are in units of 1U (0.01 m2),
# so realArea computed below is a multiplier to get the real area of a constant-alpha
# face component
def getSOCiRealAreas(SolarCells=False):
    ### accounts for both Area1U = 0.01 m2 and for fractional coverage
    ## from Boone Tate
    # Sides(2): 64% Solar Panel, 19% Aluminum Panels(outside), 17% aluminum frame rails
    # Sides(2): 57% Solar Panel, 26% aluminum panel(outside), 17% aluminum frame rails
    # Top: 60% Solar panel, 40% Aluminum frame
    # Bottom: 38% aluminum frame, 62% PCB
    ## from provided tables
    alphaSolPan = 0.92
    epsSolPan = 0.85
    CSolPan = 324
    alphaAlPanOut = 0.87
    epsAlPanOut = 0.81
    CAlPanOut = 801
    alphaAlFrameR = 0.86
    epsAlFrameR = 0.86
    CAlFrameR = 768
    alphaAlFrame = 0.08
    epsAlFrame = 0.15
    CAlFrame = 768
    alphaPCB = 0.81
    epsPCB = 0.90
    CPCB = 1544

    # should we take Kapton instead? per Table 2 in transients.pdf:
    if False:
        alphaAlPanOut = 0.11
        epsAlPanOut = 0.33

    ## nominal area for 1U side in m2
    area1U = 0.01
    realAreaFac = {}
    realAlpha = {}
    realEps = {}
    if SolarCells:
        realAreaFac['Xp'] = area1U * 0.64
        realAreaFac['Xm'] = area1U * 0.64
        realAreaFac['Yp'] = area1U * 0.57
        realAreaFac['Ym'] = area1U * 0.57
        realAreaFac['Zp'] = area1U * 0.60
        realAreaFac['Zm'] = area1U * 0.00
        for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
            realAlpha[s] = alphaSolPan
            realEps[s] = epsSolPan
    else:
        realAreaFac['Xp'] = area1U * 0.36
        realAreaFac['Xm'] = realAreaFac['Xp']
        realAreaFac['Yp'] = area1U * 0.43
        realAreaFac['Ym'] = realAreaFac['Yp']
        realAreaFac['Zp'] = area1U * 0.40
        realAreaFac['Zm'] = area1U * 1.00
        realAlpha['Xp'] = (0.19*alphaAlPanOut+0.17*alphaAlFrameR)/(0.19+0.17)
        realAlpha['Xm'] = realAlpha['Xp']
        realAlpha['Yp'] = (0.26*alphaAlPanOut+0.17*alphaAlFrameR)/(0.26+0.17)
        realAlpha['Ym'] = realAlpha['Yp']
        realAlpha['Zp'] = alphaAlFrame
        realAlpha['Zm'] = (0.62*alphaPCB+0.38*alphaAlFrame)/(0.62+0.38)
        realEps['Xp'] = (0.19*epsAlPanOut+0.17*epsAlFrameR)/(0.19+0.17)
        realEps['Xm'] = realAlpha['Xp']
        realEps['Yp'] = (0.26*epsAlPanOut+0.17*epsAlFrameR)/(0.26+0.17)
        realEps['Ym'] = realAlpha['Yp']
        realEps['Zp'] = epsAlFrame
        realEps['Zm'] = (0.62*epsPCB+0.38*epsAlFrame)/(0.62+0.38)

    return realAreaFac, realAlpha, realEps

# it works both for solar cells and other areas (called twice)
def getQsunOrbit(p, area, realAreaFac, realAlpha, eclipseFactor):
    Qsun = {}
    for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
        Qsun[s] = realAreaFac[s] * area[s] * realAlpha[s] * p['Fsun'] * eclipseFactor
    return Qsun

def getQrefOrbit(p, area, realAreaFac, realAlpha, eclipseFactor, albedoFactor):
    Qref = {}
    for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
        Qref[s] = realAreaFac[s]*area[s]*realAlpha[s]*p['fE']*albedoFactor*p['rhoE']*p['Fsun']*eclipseFactor
    return Qref

def getQIROrbit(p, area, realAreaFac, realAlpha):
    QIR = {}
    totA = 0
    totAalph = 0
    for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
        QIR[s] = realAreaFac[s] * area[s] * realAlpha[s] * p['fE'] * p['FIR']
        totA += np.mean(realAreaFac[s]) * np.mean(area[s])
        totAalph += np.mean(realAreaFac[s]) *  np.mean(area[s]) *  np.mean(realAlpha[s])
    return QIR

def getAllHeatQ(p):
    return getQsun(p), getQref(p), getQIR(p), getQdissip(p)

def sumQsForAllSides(QnotCell, Qcell, QrefNotCell, QrefCell, QIR, etaCell):
    # sum input Q for all sides (including inefficiencies of solar cells)
    QheatSol = 0*QnotCell['Xm']
    QheatRef = 0*QnotCell['Xm']
    QheatIR = 0*QnotCell['Xm']
    Qcharge = 0*QnotCell['Xm']
    for s in ['Xp', 'Xm', 'Yp', 'Ym','Zp', 'Zm']:
        QheatSol += QnotCell[s] + (1-etaCell)*Qcell[s]
        QheatRef += QrefNotCell[s]  + (1-etaCell)*QrefCell[s]
        QheatIR += QIR[s]
        Qcharge += etaCell*(Qcell[s] + QrefCell[s])
    return QheatSol, QheatRef, QheatIR, Qcharge

### 3) numerical solution from computed heating flux along the orbit

def getTsolution(alpha, Qheat, params, Tc, Pc, Niter=25):
    # translate orbital angle to time in seconds
    time = alpha / 360 * params['PorbMin'] * 60
    Tstart = 273
    for k in range(0,Niter):
        Temp, QheatTot = solveTforOrbit(Tstart, time, Qheat, params, Tc, Pc)
        Tstart = 0.5*(Tstart+Temp[-1])
    print('after', Niter, 'iterations, T tolerance = ', Tstart-Temp[-1], np.min(Temp), np.max(Temp))
    return time/60, Temp, QheatTot

def solveTforOrbit(T0, time, Qheat, params, Tc, Pc):
    T = T0 + 0*time
    emissivity = params['epsT']
    Area = params['Atot']
    C = params['C']
    mass = params['mass']
    physConst = getPhysConstants()
    sigmaSB = physConst['sigSB']
    QheatTot = 1.0*Qheat
    # direct integration with very tiny time steps assumed
    # time step controled by the orbital angle grid
    for i in range(1,np.size(time)):
        if (T[i-1] < Tc):
            QheatTot[i] = Qheat[i] + Pc
        else:
            QheatTot[i] = Qheat[i]
        tStep = time[i]-time[i-1]
        T[i] = nextTempOrbit(T[i-1], tStep, QheatTot[i], C, mass, sigmaSB, emissivity, Area)
    return T, QheatTot

def nextTempOrbit(T, t_step, QheatingTotal, c, mass, sigma, emissivity, Area):
    dTdt = (1/(c*mass))*(QheatingTotal-(sigma*(T**4)*emissivity*Area))
    return T + dTdt*t_step

### 4) new plots

def TempsPlotOrbit(time1, temp1, c1, time2="", temp2="", c2="", time3="", temp3="", c3="", outfile="", title=""):
    physConst = getPhysConstants()
    KelToCel = physConst['KelToCel']

    # Plot the Temperature wrt time
    battTmin =  0.0
    battTmax = 40.0
    # operating range for ClydeSpace batteries
    # battTmin = -10.0
    # battTmax =  50.0

    fig, ax = plt.subplots()
    fig.subplots_adjust(wspace=0.22, left=0.16, right=0.97, bottom=0.15, top=0.94)

    ax.tick_params(axis='both', which='major', labelsize=15)
    plt.plot(time1,temp1-KelToCel,label=c1, lw=3)
    if (c2 != ""):
        plt.plot(time2,temp2-KelToCel,label=c2, lw=3)
    if (c3 != ""):
        plt.plot(time3,temp3-KelToCel,label=c3, lw=3)
    plt.plot(time1, 0*time1+battTmin, lw=1, c='black')
    plt.plot(time1, 0*time1+battTmax, lw=1, c='black')
    ax.fill_between(time1, battTmin, battTmax, alpha=0.1)
    # -10 deeper during eclipse
    plt.plot([30, 66], [battTmin-20, battTmin-20], lw=1, c='black')
    ax.fill_between([30, 66], battTmin-20, battTmin, alpha=0.1)

    plt.rc('legend',fontsize=12)
    plt.legend(loc=1)
    plt.title(title)
    plt.xlabel('Time (minute)', fontsize=18)
    plt.ylabel('Temperature (C)', fontsize=18)
    plt.ylim(-25, 65.0)

    if (outfile == ""):
        name = 'figures/TempsVsOrbitTime.png'
    else:
        name = 'figures/TempTimeOrbit' + outfile + '.png'
    # plt.show()
    plt.savefig(name)
    print('made figure', name)
    plt.close("all")
    return

### PLOTTING
def TempsPowerOrbit(time1, temp1, c1, time2="", temp2="", c2="", time3="", temp3="", c3="", outfile="", title=""):

    fig, ax = plt.subplots()
    fig.subplots_adjust(wspace=0.22, left=0.16, right=0.97, bottom=0.15, top=0.94)

    ax.tick_params(axis='both', which='major', labelsize=15)
    plt.plot(time1,temp1,label=c1, lw=3)
    if (c2 != ""):
        plt.plot(time2,temp2,label=c2, lw=3)
    if (c3 != ""):
        plt.plot(time3,temp3,label=c3, lw=3)

    plt.rc('legend',fontsize=12)
    plt.legend(loc=2)
    plt.title(title)
    plt.xlabel('Time (minute)', fontsize=18)
    plt.ylabel('Power (W)', fontsize=18)
    plt.plot([0, 30], [6, 6], lw=1, c='black')
    plt.plot([66, 96], [6, 6], lw=1, c='black')
    plt.text(0, 8.0, 'EPS: 6 W')
    plt.ylim(-1, 65)

    if (outfile == ""):
        name = 'figures/PowerVsOrbitTime.png'
    else:
        name = 'figures/PowerTimeOrbit' + outfile + '.png'
    # plt.show()
    plt.savefig(name)
    print('made figure', name)
    plt.close("all")
    return

max_stars_repo_name: ivezic/CubeSats

max_issues_repo_name: ivezic/CubeSats

max_forks_repo_name: ivezic/CubeSats


text:
from bmtk.builder import NetworkBuilder
import numpy as np
import sys
import synapses
import pandas as pd
synapses.load()
syn = synapses.syn_params_dicts()

if __name__ == '__main__':
    if __file__ != sys.argv[-1]:
        inp = sys.argv[-1]
    else:
        raise Exception("no work" + str(sys.argv[-1]))

N = int(inp)

# exc_fr = float(frs[0])
# inh_fr = float(frs[1])
# print(exc_fr, inh_fr)

#N = 10

# Initialize our network

net = NetworkBuilder("biophysical")

def lognormal(m, s):
        mean = np.log(m) - 0.5 * np.log((s/m)**2+1)
        std = np.sqrt(np.log((s/m)**2 + 1))
        return max(np.random.lognormal(mean, std, 1), 0)

# Dend Excitatory: 8998.0
# Apic Excitatory: 12380.0
# Soma Inhibitory: 2285.0

# num_inh = [231, 405, 61]#[int(lognormal(56, 7.5)) for i in range(N)]
# print(num_inh)
# inh_bounds = []
# sum = 0
# for num in num_inh:
#         sum += num
#         inh_bounds.append(sum)

# num_inh_soma = [int(lognormal(61, 9)) for i in range(N)]
# print(num_inh_soma)
# inh_bounds_soma = []
# sum = 0
# for num in num_inh_soma:
#         sum += num
#         inh_bounds_soma.append(sum)

##############################Ecitatory###############################333

# num_exc = [1590, 2785]#[int(lognormal(1000, 80)) for i in range(N)]
# print(num_exc)
# exc_bounds = []
# sum = 0
# for num in num_exc:
#         sum += num
#         exc_bounds.append(sum)

# num_exc_apic = [int(lognormal(2785, 350)) for i in range(N)]
# print(num_exc_apic)
# exc_bounds_apic = []
# sum = 0
# for num in num_exc_apic:
#         sum += num
#         exc_bounds_apic.append(sum)

# num_exc_dend = [int(lognormal(1590, 200)) for i in range(N)]
# print(num_exc_dend)
# exc_bounds_dend = []
# sum = 0
# for num in num_exc_dend:
#         sum += num
#         exc_bounds_dend.append(sum)

scale_div =10

num_dend_exc = 7186 // scale_div
num_apic_exc = 10417 // scale_div
num_soma_inh = 1908 // scale_div

exc_fr_mean = 0.1
exc_fr_std = 0.5
inh_fr = 7 #* scale_div

##################################################################################
###################################Pyr Type C#####################################

net.add_nodes(N=N, pop_name='Pyrc',
    potental='exc',
    model_type='biophysical',
    model_template='hoc:L5PCtemplate',
    morphology = None)

# net.add_nodes(N=N, pop_name='Pyrc',
#     potental='exc',
#     model_type='biophysical',
#     model_template='ctdb:Biophys1.hoc',
#     model_processing='aibs_allactive',
#     #model_processing='my_allactive',
#     dynamics_params='L5Conductances.json',
#     #dynamics_params='491766131_fit.json',
#     morphology='Rbp4-Cre_KL100_Ai14-203503.04.01.01_527109145_m.swc')
#     #morphology='Rbp4-Cre_KL100_Ai14-203503.04.01.01_527109145_m.swc')

##################################################################################
###################################External Networks##############################

#print("Internal nodes built")

num_exc = (num_apic_exc + num_dend_exc) * N

# External excitatory inputs
exc_stim = NetworkBuilder('exc_stim')
exc_stim.add_nodes(N=num_exc,
                pop_name='exc_stim',
                potential='exc',
                model_type='virtual')

num_inh = num_soma_inh * N

# External inhibitory inputs
inh_stim = NetworkBuilder('inh_stim')
inh_stim.add_nodes(N=num_inh,
                pop_name='inh_stim',
                potential='exc',
                model_type='virtual')

##################################################################################
###################################Edges##########################################

def correct_cell(source, target, num_per, start):
        #import pdb; pdb.set_trace()
        sid = source.node_id - start
        tid = target.node_id
        #import pdb; pdb.set_trace()

        lower_bound = num_per * tid

        upper_bound = lower_bound + num_per

        if sid < upper_bound and sid >= lower_bound:
                #print("connecting cell {} to {}".format(sid+start,tid))
                return 1
        else:
                return None

#Inhibitory on soma.
net.add_edges(source=inh_stim.nodes(), target=net.nodes(),
                connection_rule=correct_cell,
                connection_params={'num_per': num_soma_inh , 'start':0},
                syn_weight=1,
                delay=0.1,
                dynamics_params='INT2PN.json',
                model_template=syn['INT2PN.json']['level_of_detail'],
                distance_range=[-2000.0, 2000.0],
                target_sections=['somatic'])

#start += np.sum(inh_bounds_apic)
start = 0

# Create connections between Exc --> Pyr cells

#Excitatory on basal dendrites.
net.add_edges(source=exc_stim.nodes(), target=net.nodes(),
                connection_rule=correct_cell,
                connection_params={'num_per': num_dend_exc, 'start':start},
                syn_weight=1,
                target_sections=['dend'],
                delay=0.1,
                distance_range=[-2000.0, 2000.0],
                dynamics_params='PN2PN.json',
                model_template=syn['PN2PN.json']['level_of_detail'])

#start += np.sum(exc_bounds_dend)
start += np.sum(num_dend_exc * N)

#Excitatory on apical dendrites.
net.add_edges(source=exc_stim.nodes(), target=net.nodes(),
                connection_rule=correct_cell,
                connection_params={'num_per': num_apic_exc, 'start':start},
                syn_weight=1,
                target_sections=['apic'],
                delay=0.1,
                distance_range=[-2000.0, 2000.0],
                dynamics_params='PN2PN.json',
                model_template=syn['PN2PN.json']['level_of_detail'])

# Build and save our networks
net.build()
net.save_nodes(output_dir='network')
net.save_edges(output_dir='network')

exc_stim.build()
exc_stim.save_nodes(output_dir='network')

inh_stim.build()
inh_stim.save_nodes(output_dir='network')

from bmtk.utils.reports.spike_trains import PoissonSpikeGenerator
from bmtk.utils.reports.spike_trains.spikes_file_writers import write_csv

exc_dend_frs = []
exc_apic_frs = []

exc_means = []
exc_stds = []
exc_maxs = []

for i in range(N):
        fr_mean = np.random.uniform(exc_fr_mean - 0.015, exc_fr_mean + 0.03, 1)
        #fr_mean = np.random.uniform(exc_fr_mean + 0.5, exc_fr_mean + 0.55, 1)
        dend_frs = [min(float(lognormal(fr_mean, exc_fr_std)), fr_mean+8*exc_fr_std) for _ in range(num_dend_exc)]
        apic_frs = [min(float(lognormal(fr_mean, exc_fr_std)), fr_mean+8*exc_fr_std) for _ in range(num_apic_exc)]

        frs = dend_frs + apic_frs

        exc_dend_frs += dend_frs
        exc_apic_frs += apic_frs

        exc_means.append(np.mean(frs))
        exc_stds.append(np.std(frs))
        exc_maxs.append(np.max(frs))

exc_frs = exc_dend_frs + exc_apic_frs
fr_df = pd.DataFrame()
#fr_df['gid'] = [i + num_exc+num_inh for i in range(N)]
fr_df['gid'] = [i for i in range(N)]
fr_df['fr_mean'] = exc_means
fr_df['fr_std'] = exc_stds
fr_df['fr_max'] = exc_maxs

fr_df.to_csv('frs_temp.csv', index=False)

exc_psg = PoissonSpikeGenerator(population='exc_stim')
for i in range(num_exc):
        exc_psg.add(node_ids=[i],
                firing_rate=float(exc_frs[i]/1000),
                times=(1.0*1000, 4.0*1000))
exc_psg.to_sonata('exc_stim_spikes.h5')

inh_psg = PoissonSpikeGenerator(population='inh_stim')
inh_psg.add(node_ids=range(num_inh),
        firing_rate=inh_fr/1000,
        times=(1.0*1000, 4.0*1000))
inh_psg.to_sonata('inh_stim_spikes.h5')

from bmtk.utils.sim_setup import build_env_bionet

build_env_bionet(base_dir='./',
                network_dir='./network',
                tstop=4000.0, dt = 0.1,
                report_vars=['v', 'cai'],
                current_clamp={           # Creates a step current from 500.ms to 1500.0 ms
                     'amp': [-0.5 for i in range(N)],
                     #'std': [0.0, 0.0],
                     'delay': [2000 for i in range(N)],
                     'duration': [1000 for i in range(N)],
                     'gids':"all"
                 },
                spikes_threshold=-10,
                spikes_inputs=[('exc_stim', 'exc_stim_spikes.h5'), ('inh_stim', 'inh_stim_spikes.h5')],
                components_dir='../biophys_components',
                compile_mechanisms=True)


text:
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from scipy.spatial import distance
from skimage.feature import peak_local_max
from scipy import stats,ndimage
import os
#import seaborn as sns

airyscan = False
laptop = False
files_per_batch = 18
tresh_mov = 0.2
tresh_gat = 0.25
tresh_glut1 = 0.25
tresh_glut2 = 0.5
min_dist = 3

def get_centers(peaks):
    labels, nr_objects = ndimage.label(peaks) # get all distinct features
    label_list = list(labels[np.nonzero(labels)]) # get list of feature labels
    centers = np.asarray(ndimage.center_of_mass(peaks,labels,label_list),dtype=int)   # get center of mass for all features
    return centers

def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1, length = 60, fill = '█'):
    """
    Call in a loop to create terminal progress bar
    @params:
        iteration   - Required  : current iteration (Int)
        total       - Required  : total iterations (Int)
        prefix      - Optional  : prefix string (Str)
        suffix      - Optional  : suffix string (Str)
        decimals    - Optional  : positive number of decimals in percent complete (Int)
        length      - Optional  : character length of bar (Int)
        fill        - Optional  : bar fill character (Str)
    """
    percent = ("{0:." + str(decimals) + "f}").format(100 * (iteration / float(total)))
    filledLength = int(length * iteration // total)
    bar = fill * filledLength + '-' * (length - filledLength)
    print('\r%s |%s| %s%% %s' % (prefix, bar, percent, suffix), end = '\r')
    # Print New Line on Complete
    if iteration == total:
        print()

if laptop:
    if airyscan:
        root = r'C:\Users\Hendrik\Desktop\Studium\Master\Lab_rotations\Dresbach\analysis\Airyscan'
        res = 2024
    else:
        root = r'C:\Users\Hendrik\Desktop\Studium\Master\Lab_rotations\Dresbach\analysis\Confocal'
        res = 1024
else:
    if airyscan:
        root = r'D:\Studium\Master\Lab_rotations\Dresbach\analysis\Airyscan'
        res = 2024
    else:
        root = r'D:\Studium\Master\Lab_rotations\Dresbach\analysis\Confocal'
        res = 1024

# load mover/transporter images
print('\nLoad images...')
dirnames = os.listdir(root)
ndirs = len(dirnames)
if airyscan:
    mover_im = np.zeros((res,res,files_per_batch*3))
else:
    mover_im = np.zeros((res,res,files_per_batch*3))
trans_im = mover_im.copy()
for i in range(ndirs):
    filenames = os.listdir(os.path.join(root,dirnames[i]))
    nfiles = len(filenames)
    if i%2==0:      # load mover images
        for j in range(nfiles):
            curr_file = os.path.join(root,dirnames[i],filenames[j])
            curr_im = np.array(Image.open(curr_file))
            mover_im[:,:,int(i/2)*files_per_batch+j] = curr_im
    else:
        for k in range(nfiles):
            curr_file = os.path.join(root,dirnames[i],filenames[k])
            curr_im = np.array(Image.open(curr_file))
            trans_im[:,:,int(i/2)*files_per_batch+k] = curr_im

# calculate coordinates of local maxima
print('Calculate local maxima...')
nimgs = mover_im.shape[2]
printProgressBar(0,nimgs)
mover_coord = []
trans_coord = []
for i in range(nimgs):
    curr_peaks = peak_local_max(mover_im[:,:,i],min_distance = min_dist,threshold_rel=tresh_mov,indices=False) # get local maxima
    curr_centers = get_centers(curr_peaks)
    mover_coord.append(curr_centers)
    if i < files_per_batch:
        curr_peaks = peak_local_max(trans_im[:,:,i],min_distance = min_dist,threshold_rel=tresh_gat,indices=False) # get local maxima
        curr_centers = get_centers(curr_peaks)
        trans_coord.append(curr_centers)
    elif files_per_batch <= i and i < files_per_batch*2:
        curr_peaks = peak_local_max(trans_im[:,:,i],min_distance = min_dist,threshold_rel=tresh_glut1,indices=False) # get local maxima
        curr_centers = get_centers(curr_peaks)
        trans_coord.append(curr_centers)
    else:
        curr_peaks = peak_local_max(trans_im[:,:,i],min_distance = min_dist,threshold_rel=tresh_glut2,indices=False) # get local maxima
        curr_centers = get_centers(curr_peaks)
        trans_coord.append(curr_centers)
    printProgressBar(i+1,nimgs)

#%% minimal peak distance
#  get minimal distances between mover and transporters
dist_gat = []
dist_glut1 = []
dist_glut2 = []
print('Calculate distances...')
for i in range(nimgs):
    curr_dist = distance.cdist(trans_coord[i],mover_coord[i])
    curr_min_dist = np.zeros(curr_dist.shape[0])
    for j in range(curr_dist.shape[0]):
        curr_min_dist[j] = np.min(curr_dist[j,:])
    if i < files_per_batch:
        dist_gat.append(np.median(curr_min_dist))
        #dist_gat.append(stats.mode(curr_min_dist)[0][0])
    elif files_per_batch <= i and i < files_per_batch*2:
        dist_glut1.append(np.median(curr_min_dist))
        #dist_glut1.append(stats.mode(curr_min_dist)[0][0])
    else:
        dist_glut2.append(np.median(curr_min_dist))
        #dist_glut2.append(stats.mode(curr_min_dist)[0][0])
    printProgressBar(i+1,nimgs)

dist_gat = np.array(dist_gat)
dist_glut1 = np.array(dist_glut1)
dist_glut2 = np.array(dist_glut2)

# delete outlier
dist_glut1 = np.delete(dist_glut1,(8,17))

if airyscan:
    dist_glut1 = dist_glut1[:-2]

#plt.figure('Histograms mean')
#plt.hist(dist_gat)
#plt.hist(dist_glut1)
#plt.hist(dist_glut2)

plt.figure('Boxplots median')
plt.boxplot((dist_gat,dist_glut1,dist_glut2),notch=True,labels=('vGAT','vGluT1','vGluT2'))
ax = plt.gca()
ax.set_ylim(0,25)

#%% minimal maxima distances - cumulative distribution
#  get minimal distances between mover and transporters
dist_gat = []
dist_glut1 = []
dist_glut2 = []
dist_gat = np.array(dist_gat)
dist_glut1 = np.array(dist_glut1)
dist_glut2 = np.array(dist_glut2)
print('Calculate distance distributions...')
for i in range(len(mover_coord)):
    curr_dist = distance.cdist(trans_coord[i],mover_coord[i])
    curr_min_dist = np.zeros(curr_dist.shape[0])
    for j in range(curr_dist.shape[0]):
        curr_min_dist[j] = np.min(curr_dist[j,:])
    if i < files_per_batch:
        dist_gat = np.concatenate((dist_gat,curr_min_dist))
    elif files_per_batch <= i and i < files_per_batch*2:
        dist_glut1 = np.concatenate((dist_glut1,curr_min_dist))
    else:
        dist_glut2 = np.concatenate((dist_glut2,curr_min_dist))
    printProgressBar(i+1,nimgs)

if airyscan:
    dist_glut1 = dist_glut1[:-2]

dist_glut1 = np.delete(dist_glut1,(8,17))
'''
# combined histo
plt.figure('Combined Histograms')
plt.hist(dist_gat, bins=50,label='vGAT',histtype='step',linewidth=2, cumulative=False, density=True)
plt.hist(dist_glut1,bins=50,label='vGluT1',histtype='step',linewidth=2, cumulative=False, density=True)
plt.hist(dist_glut2,bins=50,label='vGluT2',histtype='step',linewidth=2, cumulative=False, density=True)
plt.legend()
ax = plt.gca()
ax.set_xlim(0,60)
'''
# cumulative distribution function
dist_gat_sort = np.sort(dist_gat)
dist_gat_freq = np.array(range(len(dist_gat)))/float(len(dist_gat))
dist_glut1_sort = np.sort(dist_glut1)
dist_glut1_freq = np.array(range(len(dist_glut1)))/float(len(dist_glut1))
dist_glut2_sort = np.sort(dist_glut2)
dist_glut2_freq = np.array(range(len(dist_glut2)))/float(len(dist_glut2))
plt.figure('CDF')
plt.plot(dist_gat_sort,dist_gat_freq,label='vGAT')
plt.plot(dist_glut1_sort,dist_glut1_freq,label='vGluT1')
plt.plot(dist_glut2_sort,dist_glut2_freq,label='vGluT2')
plt.legend()
#ax = plt.gca()
#ax.set_xlim(0,100)
#%% below certain distance ratio (interaction ratio)
int_thresh = 5

#  get distances between mover and transporters
int_gat = []
int_glut1 = []
int_glut2 = []
for i in range(len(mover_coord)):
    curr_dist = distance.cdist(trans_coord[i],mover_coord[i])
    curr_min_dist = np.zeros(curr_dist.shape[0])
    for j in range(curr_dist.shape[0]):
        curr_min_dist[j] = np.min(curr_dist[j,:])
    if i < files_per_batch:
        int_gat.append(len(np.where(curr_min_dist<=int_thresh)[0])/curr_min_dist.shape[0])
    elif files_per_batch <= i and i < files_per_batch*2:
        try: int_glut1.append(len(np.where(curr_min_dist<=int_thresh)[0])/curr_min_dist.shape[0])
        except: pass
    else:
        int_glut2.append(len(np.where(curr_min_dist<=int_thresh)[0])/curr_min_dist.shape[0])
    printProgressBar(i+1,nimgs)

int_gat = np.array(int_gat)
int_glut1 = np.array(int_glut1)
int_glut2 = np.array(int_glut2)

int_glut1 = np.delete(int_glut1,(8,17))

#plt.figure('Histograms')
#plt.hist(dist_gat)
#plt.hist(dist_glut1)
#plt.hist(dist_glut2)

plt.figure('Boxplots')
plt.boxplot((int_gat,int_glut1,int_glut2),notch=True,labels=('vGAT','vGluT1','vGluT2'))
ax = plt.gca()
ax.set_ylim(0,1)
'''
plt.figure('Stripplots')
sns.stripplot(x=1,y=int_gat)
data = np.concatenate((int_gat[:,np.newaxis],int_glut1[:,np.newaxis],int_glut2[:,np.newaxis]),axis=1)
labels = ("vGAT",'vGluT1','vGluT2')
width=0.2
fig, ax = plt.subplots()
for i, l in enumerate(labels):
    x = np.ones(data.shape[0])*i + (np.random.rand(data.shape[0])*width-width/2.)
    ax.scatter(x, data[:,i], s=25)
    mean = data[:,i].mean()
    ax.plot([i-width/2., i+width/2.],[mean,mean], color="k")
ax.set_xticks(range(len(labels)))
ax.set_ylim(0,1)
ax.set_xticklabels(labels)
plt.show()
#plt.figure ('Barplots')
#plt.bar((1,2,3),(np.mean(int_gat),np.mean(int_glut1),np.mean(int_glut2)))
'''
#%% Statistical tests
# Normality test
p_norm = np.zeros(3)
s,p_norm[0] = stats.normaltest(dist_gat)
s,p_norm[1] = stats.normaltest(dist_glut1)
s,p_norm[2] = stats.normaltest(dist_glut2)

# Mann-Whitney-U test (nonparametric significance)
s,p_mwu_glut1_glut2 = stats.mannwhitneyu(int_glut1,int_glut2,alternative='less')
s,p_mwu_glut1_gat = stats.mannwhitneyu(int_glut1,int_gat,alternative='less')
s,p_mwu_gat_glut2 = stats.mannwhitneyu(int_gat,int_glut2,alternative='less')

# t-test
s,p_t_glut1_glut2 = stats.ttest_ind(int_glut1,int_glut2,equal_var=True)
s,p_t_glut1_gat = stats.ttest_ind(int_glut1,int_gat,equal_var=True)
s,p_t_gat_glut2 = stats.ttest_ind(int_gat,int_glut2,equal_var=True)

# KS test (distribution comparison)
D,p_dist_glut1_glut2 = stats.ks_2samp(dist_glut1,dist_glut2)
D,p_dist_glut1_gat = stats.ks_2samp(dist_glut1,dist_gat)
D,p_dist_gat_glut2 = stats.ks_2samp(dist_gat,dist_glut2)

#%% plot example images with local maxima overlaid
im_number = 5   # vGAT: 0-17, vGluT1: 18-35, vGluT2: 36-54

plt.figure('Raw Mover image min dist 3 0.2 ohne')
plt.imshow(mover_im[:,:,im_number],cmap = 'gray')
plt.plot(mover_coord[im_number][:,1],mover_coord[im_number][:,0],'r.')

plt.figure('Raw Transporter image min dist 3 0.5 ohne')
plt.imshow(trans_im[:,:,im_number],cmap = 'gray')
plt.plot(trans_coord[im_number][:,1],trans_coord[im_number][:,0],'r.')

#%% plot all images
im_count=0
fig_gat,ax_gat = plt.subplots(3,6,sharex=True,sharey=True)
for row in range(3):
    for column in range(6):
        ax_gat[row,column].imshow(trans_im[:,:,im_count],cmap='gray')
        ax_gat[row,column].plot(trans_coord[im_count][:,1],trans_coord[im_count][:,0],'r.',markersize=0.5)
        im_count += 1

im_count=0
fig_glut1,ax_glut1 = plt.subplots(3,6,sharex=True,sharey=True)
for row in range(3):
    for column in range(6):
        ax_glut1[row,column].imshow(trans_im[:,:,im_count+18],cmap='gray')
        ax_glut1[row,column].plot(trans_coord[im_count+18][:,1],trans_coord[im_count+18][:,0],'r.',markersize=0.5)
        im_count += 1

im_count=0
fig_glut2,ax_glut2 = plt.subplots(3,6,sharex=True,sharey=True)
for row in range(3):
    for column in range(6):
        ax_glut2[row,column].imshow(trans_im[:,:,im_count+36],cmap='gray')
        ax_glut2[row,column].plot(trans_coord[im_count+36][:,1],trans_coord[im_count+36][:,0],'r.',markersize=0.5)
        im_count += 1

#%% plot min_dist distribution for one image
image_num = 5   # vGAT: 0-17, vGluT1: 18-35, vGluT2: 36-54

curr_dist = distance.cdist(mover_coord[image_num],trans_coord[image_num])
gat_min_dist = np.zeros(curr_dist.shape[0])
for j in range(curr_dist.shape[0]):
        gat_min_dist[j] = np.min(curr_dist[j,:])
curr_dist = distance.cdist(mover_coord[image_num+files_per_batch],trans_coord[image_num+files_per_batch])
glut1_min_dist = np.zeros(curr_dist.shape[0])
for j in range(curr_dist.shape[0]):
        glut1_min_dist[j] = np.min(curr_dist[j,:])
curr_dist = distance.cdist(mover_coord[image_num+files_per_batch*2],trans_coord[image_num+files_per_batch*2])
glut2_min_dist = np.zeros(curr_dist.shape[0])
for j in range(curr_dist.shape[0]):
        glut2_min_dist[j] = np.min(curr_dist[j,:])

# combined histo
plt.figure('Histograms Mover-Transporter 2 cum')
plt.hist(gat_min_dist, bins=100,label='vGAT',histtype='step',linewidth=2, cumulative=True)
plt.hist(glut1_min_dist,bins=100,label='vGluT1',histtype='step',linewidth=2, cumulative=True)
plt.hist(glut2_min_dist,bins=100,label='vGluT2',histtype='step',linewidth=2, cumulative=True)
plt.legend()
ax = plt.gca()
ax.set_xlim(0,200)

#single histos
plt.figure('Histograms Mover-vGAT')
plt.hist(gat_min_dist, bins=50,label='vGAT')
ax = plt.gca()
ax.set_ylim(0,900)
plt.figure('Histogram Mover-vGluT1')
plt.hist(glut1_min_dist,bins=50,label='vGluT1')
ax = plt.gca()
ax.set_ylim(0,900)
plt.figure('Histogram Mover-vGluT2')
plt.hist(glut2_min_dist,bins=50,label='vGluT2')
ax = plt.gca()
ax.set_ylim(0,900)

reason: import numpy,from scipy


text:
#!/usr/bin/env python3
import numpy as np
from sys import exit

def isScalar(value):
    return isinstance(value, float) or isinstance(value, int)

class InputParameters:
    #########################################################################
    # Description:
    #   Input parameter handler class. All quantities defined in an input
    #   deck. Inputs checked with checkInputs function after the problem is
    #   instanciated.
    #########################################################################
    def __init__(self):
        # Options
        self.mode       = 'radhydro'
        self.limiter    = 'double minmod'

        # Geometry Specifications
        self.geometry   = 'slab'
        self.rL         = None
        self.rR         = None
        self.N          = None

        # Physical Constants
        self.gamma   = 5./3                 # thermodynamic constant for EOS
        self.C_v     = 1.                   # specific heat (jerks / cc-eV)
        self.kappa   = [577.35, 0, 1, 0]    # constants for opacity function of temp
        self.kappa_s = 0                    # scattering opacity (g / cm2)
        self.a       = 0.01372              # radiation constant (jerks / cc-keV4)
        self.c       = 299.792              # speed of light (cm )

        # Boundary Conditions
        self.bc_L_hydro_type    = None  # left boundary condition type
        self.bc_R_hydro_type    = None  # right boundary condition type
        self.bc_L_rad_type      = None  # left boundary condition type
        self.bc_R_rad_type      = None  # right boundary condition type
        self.bc_L_rad_val       = None  # left boundary value
        self.bc_R_rad_val       = None  # right boundary value

        # Initial Conditions
        self.rho    = None
        self.u      = None
        self.e      = None
        self.Er     = None

        # Iteration Parameters
        self.cfl         = None
        self.T_start     = None
        self.T_final     = None
        self.relErFactor  = None
        self.maxTimeStep = None
        self.dt          = None

    def checkInputs(self):
        # Options checks
        if (self.mode not in ['radhydro', 'hydro', 'rad']):
            exit("Need to be in radhydro, hydro, or rad mode.")
        if (self.limiter not in ['double minmod', 'minmod', 'vanLeer', 'none']):
            exit("Limiter not available. Choose 'double minmod', 'minmod', 'vanLeer', or 'none'.")
        # Geometry Checks
        if (self.geometry not in ['slab', 'cylinder', 'sphere']):
            exit("Geometry type {} not supported.".format(self.geometry))

        # Physical Constant Checks
        if (not isScalar(self.gamma) or self.gamma <= 1.0):
            exit("Need to specify gamma as a scalar larger than 1.0")

        # Initial Conditions
        if not callable(self.rho):
            exit("Need to specify rho initial condition as a function.")
        if not callable(self.u):
            exit("Need to specify u initial condition as a function.")
        if not callable(self.T):
            exit("Need to specify t initial condition as a function.")
        if not callable(self.Er):
            exit("Need to specify Er initial condition as a function.")

        # Boundary Conditions
        if self.bc_L_hydro_type not in ['fixed', 'transmissive', 'periodic']:
            exit("Need to specify fixed, transmissive, or periodic for left hydro bdry cond.")
        if self.bc_R_hydro_type not in ['fixed', 'transmissive', 'periodic']:
            exit("Need to specify fixed, transmissive, or periodic for right hydro bdry cond.")
        if self.bc_L_rad_type not in ["reflective", "source"]:
            exit("Need to specify reflective or source for left rad bdry cond.")
        if self.bc_R_rad_type not in ["reflective", "source"]:
            exit("Need to specify reflective or source for right rad bdry cond.")

        # Iteration parameters
        if (not isScalar(self.cfl) or self.cfl <= 0):
            exit("Need to specify a positive, scalar CFL factor.")
        if (not isScalar(self.T_start) or self.T_start < 0):
            exit("Need to specify a positive, scalar start time.")
        if (not isScalar(self.T_final) or self.T_final <= self.T_start):
            exit("Need to specify a scalar end time that is later than the start time.")
        if (not isScalar(self.maxTimeStep) or self.maxTimeStep > self.T_final):
            exit("Need to specify a scalar, positive maximum time step that is smaller than the final time.")
        if (not isScalar(self.dt) and self.dt != None):
            exit("Need to specify a fixed time step, or no time step.")