BlueSky fixes (OpenSky, Performance model, wypt speed, MVP)

bluesky/tools/__init__.py

@@ -1,6 +1,6 @@
 from bluesky import settings
 # Register settings defaults
-settings.set_variable_defaults(prefer_compiled=False)
+settings.set_variable_defaults(prefer_compiled=True)
 if settings.prefer_compiled:
     try:
         from . import cgeo as geo

bluesky/tools/src_cpp/README.txt.txt

@@ -0,0 +1,8 @@
+Note, to make use of this you currently still have to manually compile
+the extension by running:
+python setup.py build_ext --inplace
+in bluesky/tools/src_cpp
+and copying/moving the resulting cgeo.[so/dll] to bluesky/tools.
+
+In bluesky/tools/__init__, change the variable prefer_compiled 
+to True if not read correctly from settings.cfg

bluesky/traffic/aporasas.py

@@ -23,21 +23,11 @@ def create(self, n=1):
 
     def update(self):
         #--------- Input to Autopilot settings to follow: destination or ASAS ----------
-        # Convert the ASAS commanded speed from ground speed to TAS
-        if bs.traf.wind.winddim > 0:
-            vwn, vwe     = bs.traf.wind.getdata(bs.traf.lat, bs.traf.lon, bs.traf.alt)
-            asastasnorth = bs.traf.cr.tas * np.cos(np.radians(bs.traf.cr.trk)) - vwn
-            asastaseast  = bs.traf.cr.tas * np.sin(np.radians(bs.traf.cr.trk)) - vwe
-            asastas      = np.sqrt(asastasnorth**2 + asastaseast**2)
-        # no wind, then ground speed = TAS
-        else:
-            asastas = bs.traf.cr.tas # TAS [m/s]
-
         # Select asas if there is a conflict AND resolution is on
         # Determine desired states per channel whether to use value from ASAS or AP.
         # bs.traf.cr.active may be used as well, will set all of these channels
         self.trk = np.where(bs.traf.cr.hdgactive, bs.traf.cr.trk, bs.traf.ap.trk)
-        self.tas = np.where(bs.traf.cr.tasactive, asastas, bs.traf.ap.tas)
+        self.tas = np.where(bs.traf.cr.tasactive, bs.traf.cr.tas, bs.traf.ap.tas)
         self.alt = np.where(bs.traf.cr.altactive, bs.traf.cr.alt, bs.traf.ap.alt)
         self.vs  = np.where(bs.traf.cr.vsactive, bs.traf.cr.vs, bs.traf.ap.vs)
 

bluesky/traffic/asas/mvp.py

@@ -1,6 +1,7 @@
 ''' Conflict resolution based on the Modified Voltage Potential algorithm. '''
 import numpy as np
-from bluesky import stack
+import bluesky as bs
+from bluesky import stack, settings
 from bluesky.traffic.asas import ConflictResolution
 
 
@@ -166,19 +167,34 @@ def resolve(self, conf, ownship, intruder):
 
         # Initialize an array to store time needed to resolve vertically
         timesolveV = np.ones(ownship.ntraf) * 1e9
+
+        # Initialize an array to store if vs action is required irrespective of resolution domain
+        swvsact = np.zeros(ownship.ntraf, dtype=bool) 
 
         # Call MVP function to resolve conflicts-----------------------------------
         for ((ac1, ac2), qdr, dist, tcpa, tLOS) in zip(conf.confpairs, conf.qdr, conf.dist, conf.tcpa, conf.tLOS):
             idx1 = ownship.id.index(ac1)
             idx2 = intruder.id.index(ac2)
-
+
+            # Determine largest RPZ and HPZ of the conflict pair
+            rpz_m = np.max(conf.rpz[[idx1, idx2]] * self.resofach)
+            hpz_m = np.max(conf.hpz[[idx1, idx2]] * self.resofacv)
+
             # If A/C indexes are found, then apply MVP on this conflict pair
             # Because ADSB is ON, this is done for each aircraft separately
             if idx1 >-1 and idx2 > -1:
-                dv_mvp, tsolV = self.MVP(ownship, intruder, conf, qdr, dist, tcpa, tLOS, idx1, idx2)
+                dv_mvp, tsolV = self.MVP(ownship, intruder, conf, qdr, dist, tcpa, tLOS, idx1, idx2, rpz_m, hpz_m)
                 if tsolV < timesolveV[idx1]:
                     timesolveV[idx1] = tsolV
 
+                # Check if AC pair is in rpz and predicted to be in hpz next step
+                # hpz_m scaled with at least 1.2 resofacv for fast CL/DE traffic
+                hor_int = dist < rpz_m
+                ver_int = abs((ownship.alt[idx1] + ownship.vs[idx1]*settings.asas_dt) - \
+                              (intruder.alt[idx2] + intruder.vs[idx2]*settings.asas_dt)) < \
+                          hpz_m / self.resofacv * max(self.resofacv, 1.2)
+                swvsact[idx1] = np.logical_and(hor_int, ver_int)
+
                 # Use priority rules if activated
                 if self.swprio:
                     dv[idx1], _ = self.applyprio(dv_mvp, dv[idx1], dv[idx2], ownship.vs[idx1], intruder.vs[idx2])
@@ -235,11 +251,16 @@ def resolve(self, conf, ownship, intruder):
 
         # Determine ASAS module commands for all aircraft--------------------------
 
-        # Cap the velocity
-        newgscapped = np.maximum(ownship.perf.vmin,np.minimum(ownship.perf.vmax,newgs))
-
-        # Cap the vertical speed
-        vscapped = np.maximum(ownship.perf.vsmin,np.minimum(ownship.perf.vsmax,newvs))
+        # Compute the true airspeed from newgs (GS=TAS if no wind)
+        vwn, vwe     = bs.traf.wind.getdata(bs.traf.lat, bs.traf.lon, bs.traf.alt)
+        newtasnorth = newgs * np.cos(np.radians(newtrack)) - vwn
+        newtaseast  = newgs * np.sin(np.radians(newtrack)) - vwe
+        newtas      = np.sqrt(newtasnorth**2 + newtaseast**2)
+
+        # Cap the velocity and vertical speed
+        tascapped, vscapped, altcapped = bs.traf.perf.limits(newtas, newvs, ownship.alt, ownship.ax)
+        signvs = np.sign(newvs)
+        vscapped = np.where(np.logical_or(signvs == 0, signvs == np.sign(vscapped)), vscapped, -vscapped)
 
         # Calculate if Autopilot selected altitude should be followed. This avoids ASAS from
         # climbing or descending longer than it needs to if the autopilot leveloff
@@ -261,15 +282,17 @@ def resolve(self, conf, ownship, intruder):
         # be equal to auto pilot alt (aalt). This is to prevent a new asasalt being computed
         # using the auto pilot vertical speed (ownship.avs) using the code in line 106 (asasalttemp) when only
         # horizontal resolutions are allowed.
-        alt = alt * (1 - self.swresohoriz) + ownship.selalt * self.swresohoriz
-        return newtrack, newgscapped, vscapped, alt
+        # Additionally remaining at current alt (ownship.alt) is forced when AC are 
+        # in each others rpz and predicted to enter hpz based on vs and dt ASAS
+        alt = (alt * (1 - self.swresohoriz) + ownship.selalt * self.swresohoriz) * (1 - swvsact) \
+              + ownship.alt * swvsact
+
+        return newtrack, tascapped, vscapped, alt
 
-    def MVP(self, ownship, intruder, conf, qdr, dist, tcpa, tLOS, idx1, idx2):
+    def MVP(self, ownship, intruder, conf, qdr, dist, tcpa, tLOS, idx1, idx2, rpz_m, hpz_m):
         """Modified Voltage Potential (MVP) resolution method"""
         # Preliminary calculations-------------------------------------------------
-        # Determine largest RPZ and HPZ of the conflict pair, use lookahead of ownship
-        rpz_m = np.max(conf.rpz[[idx1, idx2]] * self.resofach)
-        hpz_m = np.max(conf.hpz[[idx1, idx2]] * self.resofacv)
+        # Use lookahead of ownship
         dtlook = conf.dtlookahead[idx1]
         # Convert qdr from degrees to radians
         qdr = np.radians(qdr)

bluesky/traffic/asas/src_cpp/README.txt.txt

@@ -0,0 +1,5 @@
+Note, to make use of this you currently still have to manually compile
+the extension by running:
+python setup.py build_ext --inplace
+in bluesky/traffic/asas/src_cpp
+and copying/moving the resulting casas.[so/dll] to bluesky/traffic/asas.

bluesky/traffic/performance/bada/perfbada.py

@@ -1,7 +1,7 @@
 """ BlueSky aircraft performance calculations using BADA 3.xx."""
 import numpy as np
 import bluesky as bs
-from bluesky.tools.aero import kts, ft, g0, vtas2cas, vcas2tas
+from bluesky.tools.aero import kts, ft, g0, vtas2cas, vcas2tas, vtas2mach
 from bluesky.traffic.performance.perfbase import PerfBase
 from bluesky.traffic.performance.legacy.performance import esf, phases, calclimits, PHASE
 from bluesky import settings
@@ -585,7 +585,7 @@ def limits(self, intent_v, intent_vs, intent_h, ax):
             self.limalt_flag, self.limvs, self.limvs_flag = calclimits(
                 vtas2cas(intent_v, bs.traf.alt), bs.traf.gs,
                 self.vmto, self.vmin, self.vmo, self.mmo,
-                bs.traf.M, bs.traf.alt, self.hmaxact,
+                vtas2mach(intent_v, bs.traf.alt), bs.traf.alt, self.hmaxact,
                 intent_h, intent_vs, self.maxthr,
                 self.thrust, self.D, bs.traf.tas,
                 self.mass, self.ESF, self.phase)

bluesky/traffic/performance/perfbase.py

@@ -73,6 +73,14 @@ def engchange(self, acid: "acid", engine_id: "txt" = ""):
             "The currently selected performance model doesn't support engine changes.",
         )
 
+    @command(name='MASS')
+    def selacmass(self, acid: 'acid', mass: 'float'):
+        """ MASS, acid, mass
+
+        Select aircraft mass for performance model."""
+        self.mass[acid] = mass
+        return True
+
     @command(name="PERFSTATS", aliases=("PERFINFO", "PERFDATA"))
     def show_performance(self, acid: "acid"):
         """Show aircraft perfromance parameters for aircraft 'acid'"""

bluesky/traffic/route.py

@@ -895,14 +895,8 @@ def direct(acidx: 'acid', wpname: 'wpinroute'):
             else:
                 alt = acrte.wpalt[wpidx]
 
-            # Check for valid Mach or CAS
-            if acrte.wpspd[wpidx] <2.0:
-                cas = mach2cas(acrte.wpspd[wpidx], alt)
-            else:
-                cas = acrte.wpspd[wpidx]
-
             # Save it for next leg
-            bs.traf.actwp.nextspd[acidx] = cas
+            bs.traf.actwp.nextspd[acidx] = acrte.wpspd[wpidx]
 
         # No speed specified for next leg
         else:

plugins/opensky.py

@@ -150,7 +150,7 @@ def update(self):
         if n_new:
             actype = [actypes.get(str(i), 'B744') for i in icao24[newac]]
             for j in range(n_new):
-                traf.cre(n_new, acid[newac][j], actype[j], lat[newac][j], lon[newac][j],\
+                traf.cre(acid[newac][j], actype[j], lat[newac][j], lon[newac][j],\
                          hdg[newac][j], alt[newac][j], spd[newac][j])
             self.my_ac[-n_new:] = True
 

plugins/windgfs.py

@@ -159,7 +159,6 @@ def loadwind(self, lat0: 'lat', lon0: 'lon', lat1: 'lat', lon1: 'lon',
             pred = 3
         elif self.hour == 24:
             ymd0 = "%04d%02d%02d" % (self.year, self.month, self.day)
-            print(ymd0)
             ymd1 = (datetime.datetime.strptime(ymd0, '%Y%m%d') + 
                     datetime.timedelta(days=1))
             self.year  = ymd1.year