% Space Based Radar
%
% A space based radar in LEO orbit around the earth is aquiring
% targets and observing them over a specified observation angle.
% After observing the one target the SBR is aquiring and observing
% another target. The slew angles of the satellite, the integrated
% observation angle and the angle from the targets to the SBR are
% computed and displayed
 
clear;
clearGraphic;
removeAll;

% loadreplaceAll("./SBR.data")

% Start of simulation
setSimEpoch(2000,3,21,19,55,0);
addEpoch('equinox',2000,3,21,12,0,0)
addEpoch('Start',2000,3,21,12,0,0)
addSatellite('SBR',6953.137,.03,63.435,0,45,0,'equinox');
addStation('Denver',I_earth,0,39.75,-104.99,1.61);
addStation('Little Rock',I_earth,0,34.74,-92.28,.087);
addStation('Chicago',I_earth,0,41.87,-87.64,.181);
 
 
%Limits:
%
% Constraint on satellite: 45 < AZ < 135 || -45 > AZ > -135 deg
SAZmin1 = 45*d2r;
SAZmax1 = 135*d2r;
SAZmin2 =-135*d2r;
SAZmax2 = -45*d2r;
 
% Constraint on target: 20 <= EL <= 70 deg
GS_ELmin = 20*d2r;
GS_ELmax = 70*d2r;
 
% Pickup Constraint: 45 <= AZ <= 85 || -45 >= AZ >= -85 deg
PCAZmin1 = 45*d2r;
PCAZmax1 = 85*d2r;
PCAZmin2 = -85*d2r;
PCAZmax2 = -45*d2r;
 
% Integration angle of observation 45 deg
cos_int_ang = cos(45*d2r);
nStations   = NumberOfNodes(1);
nMobiles    = NumberOfNodes(2);
nSatellites = NumberOfNodes(3);
iSatellites = nStations+nMobiles+1;
nNodes      = NumberOfNodes(4);
g_earth     = GeoPosition(1)*r2d;
 
for i=1:nStations
    latitude(i)=g_earth(1,i);
    longitude(i)=g_earth(2,i);
endfor
 
% open display
%id=addMapView();
id=1;
setLineColor(id,"white");
dt  = 3;    % simulation step in sec
maxt = 800; % maximum simulation time in sec
 
t      = 0:dt:maxt;
iStart = 0;
j      = 0;
jj     = 0;

setSimStepSize(dt);

for t1 = t,
    [d,r,a,el,v] = RelPosition();
    jj = jj+1;
    EE(jj) = el(3,4)*r2d;
    AA(jj)=a(4,3)*r2d;
    %if satellite is in target constraint start observation
    if iStart < 1,
        ETarget = el(1:nStations,iSatellites);
        condTarget = 0;
        for i = 1:nStations,
            if condTarget==0,
                if ETarget(i) >= GS_ELmin && ETarget(i) <= GS_ELmax,
                    eTarget = ETarget(i);
                    aSat = a(iSatellites,i);

                    %check pickup criteria
                    if ( (aSat >= PCAZmin1 && aSat <= PCAZmax1) ||
                         (aSat >= PCAZmin2 && aSat <= PCAZmax2)),
                        condTarget = i;
                    endif
                endif
            endif
        endfor
    endif
 
    if condTarget > 0,
        eTarget = el(condTarget,iSatellites);
        aTarget = a(condTarget,iSatellites);
        if eTarget >= GS_ELmin && eTarget <= GS_ELmax,
            %if target is in satellite constraint
            aSat = a(iSatellites,condTarget);
            if ((aSat >= SAZmin1 && aSat <= SAZmax1) ||
                (aSat >= SAZmin2 && aSat <= SAZmax2))
                if iStart < 1,
                    iunitvector = unitvector(aTarget,eTarget);
                    iStart = 1;
                endif
 
                uv = unitvector(aTarget,eTarget);
                cos_sw_angle = uv'*iunitvector;
 
                %within observation task
                if cos_sw_angle >= cos_int_ang,
                    j = j+1;
                    afs(j) = aSat*r2d;
                    eft(j) = eTarget*r2d;
                    D(j)   = d(condTarget,iSatellites);
                    V(j)   = r(condTarget,iSatellites);
                    T(j)   = t1;
                    OA(j)  = acos(cos_sw_angle)*r2d;
                    g_Sat  = GeoPosition(N_satellite)*r2d;
                    lonS   = g_Sat(2);
                    latS   = g_Sat(1);
                    lonT   = longitude(condTarget);
                    latT   = latitude(condTarget);
 
                   drawLine(id,lonS,latS,lonT,latT);
 
                else
                    iStart = 0;
                    condTarget=0;
                endif
            endif
        endif
    endif
 
    stepSim;

endfor
 
 
erase
 
window('111')
title('SBR: h=575 km, Observation Angle = 45 deg')
 
window('221')
plot(T,afs,'grid','red,point=3')
ylabel("AZ from Sat [deg]")

window('222')
plot(T,eft,'grid','green,point=3')
ylabel("EL from target [deg]")
 
window('223')
plot(T,D,'grid','red,point=3')
ylabel('Slant Range [km]')
 
window('224')
plot(T,OA,'grid','red,point=3')
ylabel('Observation angle [deg]')
xlabel('time [sec]')