add a few more READMEs

Author: XuyangCaoUCSD

StarPerf_MATLAB_stage/Create_Fac.m

@@ -0,0 +1,30 @@
+function Create_Fac(conid)
+disp('settings of Fac');
+global No_fac Lat Long;
+% San Diego, New York, London, Santiago, sydney, ground station 1, ground station 2,...
+Lat = [32.88 40.74 51.51 -33.45, -33.89];
+Long = [-117.23 -74.10 -0.11 -70.65, 151.20];
+
+% add Starlink ground stations
+ground_stations = readtable('../../Starlink_ground_stations.xlsx');
+num_ground_stations = height(ground_stations)
+
+for j=1:num_ground_stations
+    Lat = [Lat ground_stations{j, 1}];
+    Long = [Long ground_stations{j, 2}];
+end
+    
+No_fac=length(Long)
+for i=7:No_fac
+    info_facility=strcat('Fac',num2str(i));
+    stkNewObj('*/','Facility',info_facility);
+    lat=Lat(i);
+    long=Long(i);
+    info_facility=strcat('Scenario/Matlab_Basic/Facility/',info_facility);
+    stkSetFacPosLLA(info_facility, [lat*pi/180; long*pi/180; 0]);
+    stkConnect(conid,'SetConstraint',info_facility,'ElevationAngle Min 20');
+    num_fac(i)=i;
+    
+end
+save('Num_fac.mat','num_fac');
+end

StarPerf_MATLAB_stage/Create_LEO.m

@@ -0,0 +1,50 @@
+function [parameter] = Create_LEO(conid,path)
+% Create LEOs in STK
+% input:
+%   conid: used to connect to STK
+%   path: configuration file path of mega-constellations
+
+    global No_leo cycle No_snap tStop constellation dT;
+    parameter = readtable(path);
+    parameter = parameter.Value;
+    constellation = parameter{1,1};
+    Altitude = str2num(parameter{2,1});
+    cycle = str2num(parameter{3,1});
+    No_snap = floor(cycle/dT)+1;
+    tStop = cycle;
+    dtr = pi/180;
+    % orbit inclination
+    inc = str2num(parameter{4,1})*dtr;
+    % phase shift
+    F = str2num(parameter{5,1});
+    % number of orbits
+    leo_plane = str2num(parameter{6,1});
+    % number of satellites per orbit
+    no = str2num(parameter{7,1});
+    if (str2num(parameter{4,1}) > 80) && (str2num(parameter{4,1}) < 100)
+        raan=[0:180/leo_plane:180/leo_plane*(leo_plane-1)];
+    else
+        raan=[0:360/leo_plane:360/leo_plane*(leo_plane-1)];
+    end
+    meanAnomaly1 = [0:360/no:360/no*(no-1)];
+    raan = raan.*dtr;
+    No_leo = leo_plane*no;
+    disp('LEO:');
+    disp(No_leo);
+    for i =1:leo_plane
+        for j=1:no
+            ra = raan(i);
+            ma = (mod(meanAnomaly1(j) + 360*F/(leo_plane*no)*(i-1),360))*dtr;
+            num = j+no*(i-1);%%index of satellite-ID
+            sat_no = strcat('Sat',num2str(num));
+            stkNewObj('*/','Satellite',sat_no);
+            sat_no = strcat('*/Satellite/',sat_no);
+            stkSetPropClassical(sat_no,'J4Perturbation','J2000',0.0,tStop,dT,0,6371000 + Altitude * 10^3,0.0,inc,0.0,ra,ma);
+            num_leo(num) = num;
+        end
+    end
+    save('Num_leo.mat','num_leo','leo_plane');
+    mkdir(strcat(constellation,'\\delay'))
+end
+
+

StarPerf_MATLAB_stage/Create_delay.m

@@ -0,0 +1,58 @@
+function [ delay ] = Create_delay(position_cbf,time, inc)
+% Calculate delay between LEOs and facilities
+% input:
+%   position_cbf: used to connect to STK
+%   delay: two-dimensional matrix, delay(i,j) denotes 
+    global No_fac  No_leo constellation;
+    load('Num_leo.mat')
+    load('Num_fac.mat');
+    distance = zeros(No_fac+No_leo,No_fac+No_leo);
+    delay = zeros(No_fac+No_leo,No_fac+No_leo);
+    %calculate the distance and delay between leo and others(include leo)
+    no = No_leo/leo_plane;
+    for i=1:leo_plane
+        for j=1:no
+            cur_leo = (i-1)*no+j;
+            % if the current satellite is already the last one in the orbit
+            %, its up one is the first one
+            if j ~= no
+                up_leo = (i-1)*no+j+1;
+            else
+                 up_leo = (i-1)*no+1;
+            end
+            distance(cur_leo,up_leo) = sqrt((position_cbf{cur_leo,1}(1,time) - position_cbf{up_leo,1}(1,time))^2 + (position_cbf{cur_leo,1}(2,time) - position_cbf{up_leo,1}(2,time))^2 + (position_cbf{cur_leo,1}(3,time) - position_cbf{up_leo,1}(3,time))^2);
+            distance(up_leo,cur_leo) = distance(cur_leo,up_leo);
+            % use light speed
+            delay(cur_leo,up_leo) = distance(cur_leo,up_leo) / (3 * 10^5);
+            delay(up_leo,cur_leo) = delay(cur_leo,up_leo);
+            if i ~= leo_plane
+                right_leo = i*no+j;
+                
+            else
+                if inc > 80 && inc < 100
+                    continue;
+                else
+                    right_leo = j;
+                end
+            end
+            distance(cur_leo,right_leo) = sqrt((position_cbf{cur_leo,1}(1,time) - position_cbf{right_leo,1}(1,time))^2 + (position_cbf{cur_leo,1}(2,time) - position_cbf{right_leo,1}(2,time))^2 + (position_cbf{cur_leo,1}(3,time) - position_cbf{right_leo,1}(3,time))^2);
+            distance(right_leo,cur_leo) = distance(cur_leo,right_leo);
+            delay(cur_leo,right_leo) = distance(cur_leo,right_leo) / (3 * 10^5);
+            delay(right_leo,cur_leo) = delay(cur_leo,right_leo);
+        end
+    end
+    % each satellite will always be able to talk with the ground station
+    for i = 1:No_leo  
+        for j = No_leo + 1:No_fac+No_leo
+            distance(i,j) = sqrt((position_cbf{i,1}(1,time) - position_cbf{j,1}(1,time))^2 + (position_cbf{i,1}(2,time) - position_cbf{j,1}(2,time))^2 + (position_cbf{i,1}(3,time) - position_cbf{j,1}(3,time))^2);
+            distance(j,i) = distance(i,j);
+            delay(i,j) = distance(i,j) / (3 * 10^5);
+            delay(j,i) = delay(i,j);
+        end
+    end
+    filename = [constellation '\delay\'];
+    filename = strcat(filename,num2str(time));
+    filename = strcat(filename,'.mat');
+    save(filename,'delay') 
+end
+

StarPerf_MATLAB_stage/Create_location.m

@@ -0,0 +1,33 @@
+function [position, position_cbf]=Create_location(dT)
+global No_leo  No_fac tStart tStop No_snap Lat Long;
+load('Num_leo.mat');
+load('Num_fac.mat');
+index=1;
+position = cell(No_leo + No_fac,1);
+position_cbf = cell(No_leo + No_fac,1);
+for i=1:No_leo
+    leo_info=strcat('*/Satellite/Sat',num2str(num_leo(i)));
+    [secData, secName] = stkReport(leo_info,'LLA Position',tStart, tStop, dT);
+    lat = stkFindData(secData{1}, 'Lat');
+    long = stkFindData(secData{1}, 'Lon');
+    alt = stkFindData(secData{1}, 'Alt');
+    llapos = zeros(3,No_snap);%%[lat long high]'
+    for j = 1:No_snap
+        llapos(1,j) = llapos(1,j)+ lat(j)*180/pi;%%lat
+        llapos(2,j) = llapos(2,j)+ long(j)*180/pi;%%long
+        llapos(3,j) = llapos(3,j) + alt(j);
+    end
+    position{index} = llapos;
+    position_cbf{index} = Lla2Cbf(position{index,1});
+    index=index+1;
+    
+end
+for i=1:No_fac
+    llapos = zeros(3,No_snap);
+    llapos(1,:) = llapos(1,:)+Lat(i);
+    llapos(2,:) = llapos(2,:)+Long(i);
+    position{index} = llapos;
+    position_cbf{index} = Lla2Cbf(position{index,1});
+    index=index+1;
+end
+end

StarPerf_MATLAB_stage/Lla2Cbf.m

@@ -0,0 +1,26 @@
+function position_cbf=Lla2Cbf(position)
+%load('satData.mat');
+R=6371*10^3;
+%     No=length(position);
+%     position_cbf=cell(No,1);
+%     for s=1:No
+%         tmp_pos=position{s};
+%         r=R+ tmp_pos(3,:);% Altitude in m above the earth
+%         Theta=pi/2-tmp_pos(1,:);
+%         Phi=2*pi+tmp_pos(2,:);
+%         X=(r.*sin(Theta)).*cos(Phi);
+%         Y=(r.*sin(Theta)).*sin(Phi);
+%         Z=r.*cos(Theta);
+%         position_cbf{s}=[X; Y; Z];
+%     end
+%     position_cbf{1}(:,1:100);
+%     position_cbf{2}(:,1:100);
+
+r=R+ position(3,:);% Altitude in m above the earth
+Theta=pi/2-position(1,:)*pi/180;
+Phi=2*pi+position(2,:)*pi/180;
+X=(r.*sin(Theta)).*cos(Phi);
+Y=(r.*sin(Theta)).*sin(Phi);
+Z=r.*cos(Theta);
+position_cbf=[X; Y; Z];
+end

StarPerf_MATLAB_stage/README.md

@@ -0,0 +1,15 @@
+## LEO Constellation Construction
+
+**We adopt [StarPerf](https://github.com/SpaceNetLab/StarPerf_Simulator)'s solution for this stage**. You could go through StarPerf's README to understand the details. 
+
+### Prerequisites
+
+As mentioned in StarPerf, the solution relies on 
+
+1. [MATLAB](https://www.mathworks.com/products/matlab.html)
+2. [Ansys STK with evaluation license](https://www.ansys.com/products/missions/ansys-stk)
+
+With prerequisites installed, run the MATLAB entry-point script `build_constellation.m`, and the constellation spatial data will be generated at `constellation_outputs\`.
+
+You could also directly use those public spatial data available at [StarPerf](https://www.dropbox.com/sh/ncxf84a1m9uznm2/AABwrzHdKX6ZXsEb6FV6L3foa?dl=0) and skip this stage.
+

StarPerf_MATLAB_stage/build_constellation.m

@@ -0,0 +1,50 @@
+clear all; close all; clc;
+global cycle No_snap No_fac No_leo tStart tStop dT constellation
+dT = 1.0;
+tStart = 0;
+dtr = pi/180;
+rtd = 180/pi;
+remMachine = stkDefaultHost;
+delete(get(0,'children'));
+conid=stkOpen(remMachine);
+
+scen_open = stkValidScen;
+if scen_open == 1
+    rtn = questdlg('Close the current scenario?');
+    if ~strcmp(rtn,'Yes')
+        stkClose(conid)
+    else
+        stkUnload('/*')
+    end
+end
+
+disp('Create a new scenario');
+stkNewObj('/','Scenario','Matlab_Basic');
+disp('Set scenario time period');
+stkSetTimePeriod('1 Dec 2019 0:00:00.0','1 Dec 2019 10:00:00.0','GREGUTC');
+stkSetEpoch('1 Dec 2019 0:00:00.0','GREGUTC');
+cmd1 = ['SetValues "1 Dec 2019 0:00:00.0" ' mat2str(dT)];
+cmd1 = [cmd1 ' 0.1'];
+rtn = stkConnect(conid,'Animate','Scenario/Matlab_Basic',cmd1);
+rtn = stkConnect(conid,'Animate','Scenario/Matlab_Basic','Reset');
+disp('Set up the propagator and nodes for the satellites');
+% spanning satellites
+[parameter] = Create_LEO(conid,'../constellation_params/parameter.xlsx');
+% spanning ground stations
+Create_Fac(conid);
+inc = str2num(parameter{4,1})*dtr;
+
+disp('save position info');
+[position, position_cbf]=Create_location(dT);
+filename = ['..\constellation_outputs\' constellation '\position.mat'];
+save(filename,'position','position_cbf');
+disp('save delay info');
+% for each cycle, caculate the delay matrix
+for t = 1:cycle
+    [delay] = Create_delay(position_cbf,t,inc);
+end
+stkExec( conid, 'Animate Scenario/Matlab_Basic  Reset' );
+stkExec( conid, 'Animate Scenario/Matlab_Basic  Start' );
+
+stkClose(conid)
+stkClose

constellation_params/README.md

@@ -0,0 +1,20 @@
+## LEO Constellation Parameter Specification
+
+(Adopting StarPerf's Framework)
+
+Please specify the following key parameters in `parameter.xlsx`, following the self-explanatory format:
+
+* Name: just a unique identifier of the constellation.
+* Altitude: the altitude of satellites. See the visualization.
+* Cycle: the orbital period (in seconds). I.e., the time it takes for a satellite to go back to the original spatial point.
+* Inclination: the angle between orbits and the equator. See the visualization.
+* Phase shift: the amount of shift of same-index satellites on two adjacent orbits. Usually it is 0. 
+* Number of orbits: as name suggests. See the visualization.
+* Number of satellites: the number of satellites on each orbit. Therefore, number of orbits × number of satellites = total satellite number in the constellation.
+
+Visualization of Starlink Shell 1 constellation:
+
+![Image: starlink_visualization.jpg](https://github.com/XuyangCaoUCSD/LeoEM/blob/main/constellation_params/starlink_visualization.jpg)
+
+Then we can move to the stage 1: `StarPerf_MATLAB_stage/`.
+

constellation_params/constellation.json

@@ -0,0 +1 @@
+{"constellation": [{"id": "0", "name": "starlink", "inclination": 53, "altitude": 550, "num_satellite": 1584, "phase_shift": 1, "orbit_num": 24, "satellite_num": 66, "link_band": "Ka"}, {"id": "1", "name": "starlink", "inclination": 89.7, "altitude": 1200, "num_satellite": 720, "phase_shift": 0, "orbit_num": 18, "satellite_num": 40, "link_band": "Ka"}, {"id": "2", "name": "telasat(polar)", "inclination": 99.5, "altitude": 1000, "num_satellite": 72, "phase_shift": 0, "orbit_num": 6, "satellite_num": 12, "link_band": "Ka"}, {"id": "3", "name": "telasat(inclined)", "inclination": 37.4, "altitude": 1200, "num_satellite": 50, "phase_shift": 0, "orbit_num": 5, "satellite_num": 10, "link_band": "Ka"}]}

constellation_params/parameter.xlsx

@@ -0,0 +1,47 @@
+## Route Computation
+
+We pre-compute the time-varying route between two endpoints using the time-varying spatial data calculated from the previous stage.
+
+### Usage:
+
+```bash
+route_stage % python3 precompute_path_15s_isl.py|precompute_path_15s_bp.py constellation output_filename sat_num cycle latitude_1 longitude_1 latitude_2 longitude_2 depression_angle elevation_angle
+```
+
+* To use inter-satellite laser link, choose `precompute_path_15s_isl.py`. To use bent-pipe radio link, choose `precompute_path_15s_bp.py`.
+* `constellation`: the LEO constellation name. The program will seek spatial data in `../constellation_outputs/$constellation/`. 
+* `output_filename`: the output filename. The file will be generated in `../precomputed_paths/`.
+* `sat_num`: total satellite number in the constellation. You should know it in stage 1.
+* `cycle`: the orbital period length (in seconds). You should know it in stage 1.
+* `latitude_1`: the latitude of the first endpoint you choose. In degree.
+* `longitude_1`: the longitude of the first endpoint you choose. In degree.
+* `latitude_2`: the latitude of the second endpoint you choose. In degree.
+* `longitude_2`: the longitude of the second endpoint you choose. In degree.
+* `depression_angle`: please refer to the figure in Principle subsection. In degree.
+* `elevation_angle`: please refer to the figure in Principle subsection. In degree.
+
+Also specify all the rely ground stations' locations in `../ground_stations.xlsx`, if you use bent-pipe links. The given one has Starlink's *discovered* ground stations from [here](https://satellitemap.space/).
+
+For example, to compute the dynamic ISL routes between San Diego and New York City using Starlink Shell 1 (whose spatial data has been computed in stage 1):
+
+```bash
+route_stage % python3 precompute_path_15s_isl.py Starlink Starlink_SD_NY_15_ISL_path.log 1584 5731 32.881 -117.237 40.845 -73.932 44.85 40
+```
+
+To compute the dynamic bent-pipe routes between San Diego and Seattle using Starlink Shell 1:
+
+```bash
+route_stage % python3 precompute_path_15s_bp.py Starlink Starlink_SD_SEA_15_BP_path.log 1584 5731 32.881 -117.237 47.608 -122.335 44.85 40
+```
+
+### Principle
+
+Explanation of mathematical function `getCoverageLimitL` in `utility.py`:
+
+![Image: cover.jpg](https://github.com/XuyangCaoUCSD/LeoEM/blob/main/route_stage/cover.jpg)
+
+
+![Image: elevation.jpg](https://github.com/XuyangCaoUCSD/LeoEM/blob/main/route_stage/elevation.jpg)
+
+
+