#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import print_function
import numpy as np

# Constants
c = 2.9979e8 # speed of light [m/s]
kB = 1.38064852e-23 # Boltzmann's constant [J/K]
rEarth = 6378e3 # [m]
muEarth = 3.986004418e14 # standard gravitational parameter [m^3/s^2]

# Broadband LEO system parameters
Ncells = 400e3
# Parameters derived from FCC filings
# SAT-LOA-20170726-00110
#     https://licensing.fcc.gov/myibfs/download.do?attachment_key=1380782
# SAT-MOD-20181108-00083
#     https://licensing.fcc.gov/myibfs/download.do?attachment_key=1569860
# SAT-MOD-20190830-00087
#     https://licensing.fcc.gov/myibfs/download.do?attachment_key=1877844
# Many of these filings have attached Microsoft Access database files (.mdb)
# tabulating orbit, beam, channel, and other properties of the system.
Nsats = 10e3 # Rounded up to include planned VLEO constellation from other filings
Nchannels = 76 # number of orthogonal TX Service channels
Nbeams = 15 # number of TX beams with Service channels
#
# According to SAT-MOD-20190830-00087, there are 15 TX beams per SV
# capable of providing service to users:
#
# Beam    Frequency       PFD [dB(W/m^2/MHz)]
# TX1     10700-12700     -122.0
# TX7     10700-12700     -122.0
# TX19    10700-12700     -122.0
# 
# TX2     17800-18600     -116.3
# TX8     17800-18600     -116.3
# TX14    17800-18600     -116.3
# TX20    17800-18600     -116.3
# 
# TX3     18800-19300     -116.3
# TX9     18800-19300     -116.3
# TX15    18800-19300     -116.3
# TX21    18800-19300     -116.3
# 
# TX4     19700-20200     -129.2
# TX10    19700-20200     -129.2
# TX16    19700-20200     -129.2
# TX22    19700-20200     -129.2
#
# Each of these beam is associated with a number of 50 MHz channels.  The number
# of channels and their frequencies vary from beam to beam; each beam has
# between 10 and 40 channels.  It is not entirely clear whether each beam can
# transmit on all of its channels simultaneously, or on only a subset at a time.
#
# We define a "beam-channel" to be the unit of space-frequency reservation on
# the SV.  The number of beam-channels, Nbc, is the sum over beams on an SV of
# the number of channels that this beam can transmit on simultaneously.
# 
# Based on earlier filings that listed 250 MHz downlink channels, it seems
# likely that Nbc is at least 5*Nbeams.  It may be as large as 264.
# 
Nbc = 264 # sum over TX beams of number of Service channels that can transmit at once
D = 45e3 # service cell diameter [m]
f = 12e9 # [Hz]
lambda_ = c/f # [m]
altitudeMin = 550e3 # [m]
altitudeAll = np.array([550e3, 1150e3]) # [m]
Grx = 33.2 # receive antenna gain [dBi]
# Parameters based on conjecture (reasoning in paper)
Nneighbors = 6
Tsetup = 5e-3 # [s]
Tswitch = 100e-6 # [s]
TswitchTX = TswitchRX = Tswitch
# Parameters from Supplement B.
PARmin = 9.6 # minimum peak-to-average power ratio for energy budgeting
# Parameters from ITU Radio Rules
PFD = -122.0 # receive power flux density [dB(W/m^2/MHz)]

# Friis equation
Ae  = (lambda_**2/4/np.pi) * 10**(.1 * Grx) # receive antenna effective area [m^2]
W   = 50e6 # receive bandwidth [Hz]
Prx = Ae * (W/1e6) * 10**(.1*PFD) # receive signal power [W]
Trx = 273 * 10**0.6 # receiver noise temperature [K]
Pnoise = kB * Trx * W # receive noise power [W]
SNR = Prx/Pnoise
# If we had a Shannon-capacity-achieving FEC scheme,
Cshannon = W * np.log2(1 + SNR)
# Assume instead a gap to capacity of 0.5 dB; then
Rturbo = W * np.log2(1 + SNR * 10**-(.1 * 0.5))
# System downlink capacity is then bounded by
totalDownlink = Rturbo * Nsats * Nbc

# Fused LEO system parameters
n = 5
Tburst = 500e-6 # [s]
Tperiod = 1. # [s]
minElevation = 40 # [degrees]
Tquantize = 1e-6
serviceLatitudeRange = (-60., 60.) # [degrees]
serviceLatMax = max(abs(np.array(serviceLatitudeRange)))

# Check that our round number of cells is not too small
sinServiceLatitude = [np.sin(np.deg2rad(lat)) for lat in serviceLatitudeRange]
serviceArea = 2*np.pi*rEarth**2 * (sinServiceLatitude[1] - sinServiceLatitude[0])
cellArea = 3*3**.5/8 * D**2
assert serviceArea/cellArea <= Ncells

# Compute TX, RX, mean, steering, and power reservations
TXreservationMax = Ncells * (Tburst + (n-1)*(Tburst+2*TswitchTX) * Nbc/Nbeams) / (Nbc * Nsats * Tperiod)
Tsweepmax = D/c * np.cos(np.deg2rad(minElevation))
T1 = Tburst + Tsweepmax
RXreservationMax = ((Nneighbors+1)*n*T1 + 2*(n-1)*TswitchRX) / (Nchannels * Tperiod)
downlinkReservationMax = TXreservationMax + RXreservationMax
downlinkReservationRound = downlinkReservationMax
downlinkRemaining = 1-downlinkReservationMax
setupReservationMax = ((n-1)*Ncells*Tsetup)/(Nbeams * Nsats * Tperiod)
setupReservationRound = setupReservationMax
pointingLossMaxBefore = 0.5 # dB
pointingLossMaxAfter = pointingLossMaxBefore * (1 + 2 * setupReservationMax)
pointingLossDelta = pointingLossMaxAfter - pointingLossMaxBefore
powerReservationMax = n*Ncells*Tburst*PARmin / (Nbc * Nsats * Tperiod)

# Usage of downlink bandwidth
broadbandEquivalentService = totalDownlink * downlinkReservationMax / Ncells
broadbandEquivalentServiceBytes = broadbandEquivalentService/8

# Steering rates
rOrbit = rEarth + altitudeAll
vOrbit = (muEarth/rOrbit)**.5
omegaOrbit = vOrbit / rOrbit # Rate of absolute rotation of nadir-pointing SV
omegaEarth = 2*np.pi/86400 # Rate of rotation of Earth [rad/s]
# For an SV in a polar orbit, add up east-west and north-south components in
# quadrature.  The north-south angular rate is given by the rate of north-south
# rotation of the SV-user line, less the rate of rotation of the SV (as it
# maintains nadir orientation).  The east-west angular rate is given by the rate
# of east-west rotation of the SV-user line due to Earth rotation.
omegaLineOfSight = ((vOrbit/altitudeAll - omegaOrbit)**2 + (rEarth*omegaEarth/altitudeAll)**2)**.5
broadbandSteeringRate = np.rad2deg(omegaLineOfSight.max())
# Using contours from the SAT-MOD-20190830-00087 filing, the FWHM TX beam spot
# size at boresight may be estimated as
# 
# Shell     Band    Linear size     Angular size
#  550km    KU      28.7 km         3.0°
#           KA      13.2 km         1.4°
# 1150km    KU      40.9 km         2.0°
#           KA      25.6 km         1.3°
cellLinearSizeMax = 40.9e3 # [m]
cellAreaMax = (3*3**.5/8) * cellLinearSizeMax**2

# Usage of uplink bandwidth
minElevationAngle = np.r_[25, 40] # [degrees]
thetaMin = np.deg2rad(minElevationAngle)
r,a = rEarth, altitudeAll
# Central angle subtended by arc from point directly below SV
# to point where SV elevation angle drops to minElevationAngle.
phiMax = np.arccos(r/(r+a) * np.cos(thetaMin)) - thetaMin
# Distance between SV and point where SV elevation is minElevationAngle.
serviceRangeMax = ((r+a)**2 + r**2 - 2*r*(r+a)*np.cos(phiMax))**.5
# Bound on time-of-flight
Tflightmax = serviceRangeMax.max() / c
# Using e.g. arithmetic coding, the assignment tuple may be packed into this
# number of bits:
assignmentSizeBits = int(np.ceil(np.log2(
    Nsats * Nchannels * (Tperiod / Tquantize) *
    (Tflightmax / Tquantize) *
    (Tsweepmax / Tquantize)
)))
# Users need to know where to point their antennas, and so SVs must
# know coarse ephemeris for each secondary assignment.  The amount of
# accuracy needed is "enough to point the beam to within half a degree during
# the next second or two".
# So the number of bits that must be sent is, at most, a few times the log2 of
# the number of 0.5° spot beams that fit in the sky.
spotBeamSolidAngle = 4*np.pi*np.sin(np.deg2rad(0.5))**2
skySolidAngle = 4*np.pi*np.sin(np.deg2rad(90 - minElevationAngle.min()))**2
secondaryEphemerisSizeBits = 2 * np.log2(skySolidAngle / spotBeamSolidAngle)
# Number of uplink bits needed for scheduling
assignmentUplinkCost = (2*n-1)*Ncells*assignmentSizeBits
uplinkCost = assignmentUplinkCost + (n-1)*Ncells*secondaryEphemerisSizeBits
uplinkPeriod = 15. # how often new GNSS schedules are uploaded; [s]
uplinkRatePerSV = uplinkCost/Nsats/uplinkPeriod

from supplement_common import export

export('n')
export('Ncells')
export('Nsats')
export('Nchannels')
export('Nbeams')
export('Nbc')
export('Nneighbors')
export('W', scale=1e-6, precision=0, siunit='\MHz')
export('D', scale=1e-3, siunit='\kilo\meter')
export('Tburst', scale=1e6, siunit='\micro\second')
export('Tsetup', scale=1e3, siunit='\milli\second')
export('Tswitch', scale=1e6, siunit='\micro\second')
export('Tperiod', siunit='\second')
export('minElevation', siunit='\degree')
export('Tquantize', scale=1e6, siunit='\micro\second')
export('TXreservationMax', scale=1e2, precision=2, siunit='\percent')
export('RXreservationMax', scale=1e2, precision=2, siunit='\percent')
export('downlinkReservationMax', scale=1e2, precision=2, siunit='\percent')
export('downlinkReservationRound', scale=1e2, precision=1, siunit='\percent')
export('downlinkRemaining', scale=1e2, precision=1, siunit='\percent')
export('setupReservationMax', precision=1, scale=1e2, siunit='\percent')
export('setupReservationRound', precision=0, scale=1e2, siunit='\percent')
export('pointingLossMaxAfter', precision=1, siunit='\dB')
export('pointingLossDelta', precision=1, siunit='\dB')
export('powerReservationMax', scale=1e2, precision=2, siunit='\percent')
export('serviceLatMax', siunit='\degree')
export('PARmin', precision=1, siunit='')
if broadbandEquivalentService > 1e6:
    export('broadbandEquivalentService', scale=1e-6, precision=1, siunit='Mbps')
else:
    export('broadbandEquivalentService', scale=1e-3, precision=0, siunit='kbps')
if broadbandEquivalentServiceBytes > 1e6:
    export('broadbandEquivalentServiceBytes', scale=1e-6, precision=1, siunit=r'\mega\byte\per\second')
else:
    export('broadbandEquivalentServiceBytes', scale=1e-3, precision=0, siunit=r'\kilo\byte\per\second')
export('broadbandSteeringRate', precision=2, siunit='\degree\per\second')
export('uplinkCost', precision=0, scale=2**-23, siunit=r'\mebi\byte')
export('assignmentUplinkCost', precision=0, scale=2**-23, siunit=r'\mebi\byte')
export('uplinkRatePerSV', precision=1, scale=2**-10, siunit='kbps')
export('assignmentSizeBits')
export('PFD', precision=1, siunit='\dB(\W\per\meter^2\per\MHz)')
export('secondaryEphemerisSizeBits')
export('cellAreaMax', precision=-1, scale=1e-6, siunit='\km^2')