diff --git a/wetb/flex/__init__.py b/wetb/flex/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..9bdc6ddbf9c798d100054b06b33e5d956b703bcf
--- /dev/null
+++ b/wetb/flex/__init__.py
@@ -0,0 +1,7 @@
+from ._io import load, read_sensor_info
+from wetb import gtsdf
+import numpy as np
+
+class Dataset(gtsdf.Dataset):
+ def __init__(self, filename, name_stop=8,dtype=np.float):
+ self.time, self.data, self.info = load(filename, name_stop,dtype)
\ No newline at end of file
diff --git a/wetb/flex/_io.py b/wetb/flex/_io.py
new file mode 100644
index 0000000000000000000000000000000000000000..f9419f8b409bacc87ba8c6f42adcbf8f1e631934
--- /dev/null
+++ b/wetb/flex/_io.py
@@ -0,0 +1,179 @@
+'''
+Created on 11. apr. 2017
+
+@author: mmpe
+'''
+import struct
+import numpy as np
+import os
+
+def load(filename, name_stop=8, dtype=np.float):
+ """
+ Parameters
+ ----------
+ - filename
+ - name_stop : int or str
+ if int: Number of characters for name
+ if str: name-description delimiter, e.g. " "
+ - dtype
+ """
+ if isinstance(filename, str):
+ fid = open(filename,'rb')
+ elif hasattr(filename, "name"):
+ fid = filename
+ filename = fid.name
+ try:
+ _ = struct.unpack('i', fid.read(4))
+ _ = struct.unpack('i', fid.read(4))
+ title = fid.read(60).strip()
+
+ _ = struct.unpack('i', fid.read(4))
+ _ = struct.unpack('i', fid.read(4))
+ no_sensors = struct.unpack('i', fid.read(4))[0]
+
+ sensor_numbers = [struct.unpack('i', fid.read(4))[0] for _ in range(no_sensors)]
+ _ = struct.unpack('i', fid.read(4))
+ _ = struct.unpack('i', fid.read(4))
+ time_start = struct.unpack('f', fid.read(4))[0]
+ time_step = struct.unpack('f', fid.read(4))[0]
+ scale_factors = np.array([struct.unpack('f', fid.read(4))[0] for _ in range(no_sensors)], dtype=np.double)
+ data = np.fromstring(fid.read(), 'int16').astype(dtype)
+ finally:
+ fid.close()
+
+# if title.isalnum():
+# self.set_info(title, "", self.filename)
+# else:
+# self.set_info(os.path.basename(self.name), "", self.name)
+ try:
+ data = data.reshape(len(data) // no_sensors, no_sensors)
+ except ValueError:
+ raise ValueError("The number of data values (%d) is not divisible by the number of sensors (%d)" % (len(data), no_sensors))
+
+ for i in range(data.shape[1]):
+ data[:, i] *= scale_factors[i]
+ no_scans = data.shape[0]
+
+
+ # Load sensor file if exists
+
+
+ sensor_filename = os.path.join(os.path.dirname(filename), "sensor")
+ sensor_info = {info[0]:info[1:] for info in read_sensor_info(sensor_filename, name_stop) }
+
+ # set gain and offset for "Time"
+ gains = []
+ offsets = []
+ names, units, descriptions = [], [], []
+
+ for sensor_nr in sensor_numbers:
+
+ name, unit, description, gain, offset = sensor_info.get(sensor_nr, ["Attribute %d"%sensor_nr, '-','',1,0])
+ if sensor_nr==1 and name=="Time" and unit=="s":
+ data = data[:,1:]
+ continue
+ names.append(name)
+ units.append(unit)
+ descriptions.append(description)
+ gains.append(gain)
+ offsets.append(offset)
+
+ time = np.arange(time_start, time_start + data.shape[0] * time_step, time_step).reshape((no_scans, 1))
+ #data = np.concatenate((time, data), axis=1)
+ #gains = np.r_[1,gains]
+ #offsets = np.r_[0,offsets]
+ # self[:]*=self.gains
+ # self[:]+=self.offsets
+ info = {"name": title,
+ "description": filename,
+ "attribute_names": names,
+ "attribute_units": units,
+ "attribute_descriptions": descriptions}
+ return time, data, info
+
+
+def read_sensor_info(sensor_file, name_stop=8):
+ """
+ Parameters
+ ----------
+ - sensor_file
+ - name_stop : int or str
+ if int: Number of characters for name
+ if str: name-description delimiter, e.g. " "
+ """
+
+ if hasattr(sensor_file, 'readlines'):
+ sensor_info_lines = sensor_file.readlines()[2:]
+ else:
+ with open(sensor_file, encoding="utf-8") as fid:
+ sensor_info_lines = fid.readlines()[2:]
+ sensor_info = []
+ for line in sensor_info_lines:
+ # while " " in line:
+ # line = line.replace(" "," ")
+ line = line.strip().split()
+ nr = int(line[0])
+ gain = float(line[1])
+ offset = float(line[2])
+ unit = line[5]
+ if isinstance(name_stop,int):
+ name_desc = " ".join(line[6:])
+ name = name_desc[:name_stop].strip()
+ description = name_desc[name_stop:].strip()
+ elif isinstance(name_stop,str):
+ name_desc = (" ".join(line[6:])).split(name_stop)
+ name = name_desc[0].strip()
+ description = name_stop.join(name_desc[1:]).strip()
+
+
+ sensor_info.append((nr, name, unit, description, gain, offset))
+ return sensor_info
+
+# def save(dataset, filename):
+# ds = dataset
+# # Write int data file
+# f = open(filename, 'wb')
+# time_att = ds.basis_attribute()
+# sensors = [s for s in ds.attributes() if s is not time_att]
+#
+# if isinstance(ds, FLEX4Dataset):
+# data = ds[:] # (ds[:]-ds.offsets)/ds.gains
+# else:
+# data = ds[:]
+# if time_att.index != -1: # time_att may be "Index" with index=-1 if "Time" not exists
+# data = np.delete(data, time_att.index, axis=1)
+# f.write(struct.pack('ii', 0, 0)) # 2x empty int
+# title = ("%-60s" % str(ds.name)).encode()
+# f.write(struct.pack('60s', title)) # title
+# f.write(struct.pack('ii', 0, 0)) # 2x empty int
+# ns = len(sensors)
+# f.write(struct.pack('i', ns))
+# f.write(struct.pack('i' * ns, *range(1, ns + 1))) # sensor number
+# f.write(struct.pack('ii', 0, 0)) # 2x empty int
+# time = ds.basis_attribute()
+# f.write(struct.pack('ff', time[0], time[1] - time[0])) # start time and time step
+#
+# scale_factors = np.max(np.abs(data), 0) / 32000
+# f.write(struct.pack('f' * len(scale_factors), *scale_factors))
+# # avoid dividing by zero
+# not0 = np.where(scale_factors != 0)
+# data[:, not0] /= scale_factors[not0]
+# #flatten and round
+# data = np.round(data.flatten()).astype(np.int16)
+# f.write(struct.pack('h' * len(data), *data.tolist()))
+# f.close()
+#
+# # write sensor file
+# f = open(os.path.join(os.path.dirname(filename), 'sensor'), 'w')
+# f.write("Sensor list for %s\n" % filename)
+# f.write(" No forst offset korr. c Volt Unit Navn Beskrivelse------------\n")
+# sensorlineformat = "%3s %.3f %.3f 1.00 0.00 %7s %-8s %s\n"
+#
+# if isinstance(ds, FLEX4Dataset):
+# gains = np.r_[ds.gains[1:], np.ones(ds.shape[1] - len(ds.gains))]
+# offsets = np.r_[ds.offsets[1:], np.zeros(ds.shape[1] - len(ds.offsets))]
+# sensorlines = [sensorlineformat % ((nr + 1), gain, offset, att.unit[:7], att.name.replace(" ", "_")[:8], att.description[:512]) for nr, att, gain, offset in zip(range(ns), sensors, gains, offsets)]
+# else:
+# sensorlines = [sensorlineformat % ((nr + 1), 1, 0, att.unit[:7], att.name.replace(" ", "_")[:8], att.description[:512]) for nr, att in enumerate(sensors)]
+# f.writelines(sensorlines)
+# f.close()
\ No newline at end of file
diff --git a/wetb/flex/tests/__init__.py b/wetb/flex/tests/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/flex/tests/test_files/test1/sensor b/wetb/flex/tests/test_files/test1/sensor
new file mode 100644
index 0000000000000000000000000000000000000000..88e89f63990b9e53bef2ba37ebcd5abef5083e88
--- /dev/null
+++ b/wetb/flex/tests/test_files/test1/sensor
@@ -0,0 +1,9 @@
+Sensor list for C:/mmpe/programming/python/WindEnergyToolbox/wetb/flex/tests/test_files/test1/test.int
+ No forst offset korr. c Volt Unit Navn Beskrivelse------------
+ 1 1.000 0.000 1.00 0.00 m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos 0.75, 0.00, -40.75
+ 2 1.000 0.000 1.00 0.00 m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos 1.00, 0.00, -39.00
+ 3 1.000 0.000 1.00 0.00 m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos 7.00, 0.00, -33.00
+ 4 1.000 0.000 1.00 0.00 rad/s Omega Rotor speed
+ 5 1.000 0.000 1.00 0.00 m Ae_pos_x Aero position x of blade 1 at radius 20.50, global coo.
+ 6 1.000 0.000 1.00 0.00 m Ae_pos_y Aero position y of blade 1 at radius 20.50, global coo.
+ 7 1.000 0.000 1.00 0.00 m Ae_pos_z Aero position z of blade 1 at radius 20.50, global coo.
diff --git a/wetb/flex/tests/test_files/test1/test.int b/wetb/flex/tests/test_files/test1/test.int
new file mode 100644
index 0000000000000000000000000000000000000000..5df18a5d5b28c9b45ec0b48786925a5fb7d9a306
Binary files /dev/null and b/wetb/flex/tests/test_files/test1/test.int differ
diff --git a/wetb/flex/tests/test_files/test_sensor_info/sensor b/wetb/flex/tests/test_files/test_sensor_info/sensor
new file mode 100644
index 0000000000000000000000000000000000000000..593044b9a59183b02162c7acc03652b2fa20f864
--- /dev/null
+++ b/wetb/flex/tests/test_files/test_sensor_info/sensor
@@ -0,0 +1,69 @@
+ Version ID : HAWC2MB 12.3 (beta -rev_TJ6)
+ No forst offset korr. c Volt Unit Navn Beskrivelse---------------
+ 1 1.0 0.0 0.00 1.0 s Time Time
+ 2 1.0 0.0 0.00 1.0 deg bea1 an shaft_rot angle
+ 3 1.0 0.0 0.00 1.0 rpm bea1 an shaft_rot angle speed
+ 4 1.0 0.0 0.00 1.0 deg bea2 an pitch1 angle
+ 5 1.0 0.0 0.00 1.0 deg/s bea2 an pitch1 angle speed
+ 6 1.0 0.0 0.00 1.0 deg bea2 an pitch2 angle
+ 7 1.0 0.0 0.00 1.0 deg/s bea2 an pitch2 angle speed
+ 8 1.0 0.0 0.00 1.0 deg bea2 an pitch3 angle
+ 9 1.0 0.0 0.00 1.0 deg/s bea2 an pitch3 angle speed
+ 10 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 11 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 12 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 13 1.0 0.0 0.00 1.0 kN Fx coo: Force Fx Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 14 1.0 0.0 0.00 1.0 kN Fy coo: Force Fy Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 15 1.0 0.0 0.00 1.0 kN Fz coo: Force Fz Mbdy:towertop nodenr: 1 coo: tower tower top -1: below top mass
+ 16 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:tower nodenr: 1 coo: tower tower base flange
+ 17 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:tower nodenr: 1 coo: tower tower base flange
+ 18 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:tower nodenr: 1 coo: tower tower base flange
+ 19 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:towertop nodenr: 2 coo: towertop yaw bearing
+ 20 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:towertop nodenr: 2 coo: towertop yaw bearing
+ 21 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:towertop nodenr: 2 coo: towertop yaw bearing
+ 22 1.0 0.0 0.00 1.0 kN Fx coo: Force Fx Mbdy:towertop nodenr: 2 coo: towertop yaw bering
+ 23 1.0 0.0 0.00 1.0 kN Fy coo: Force Fy Mbdy:towertop nodenr: 2 coo: towertop yaw bering
+ 24 1.0 0.0 0.00 1.0 kN Fz coo: Force Fz Mbdy:towertop nodenr: 2 coo: towertop yaw bering
+ 25 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:shaft nodenr: 4 coo: shaft main bearing
+ 26 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:shaft nodenr: 4 coo: shaft main bearing
+ 27 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:shaft nodenr: 4 coo: shaft main bearing
+ 28 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:blade1 nodenr: 3 coo: blade1 blade 1 root
+ 29 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:blade1 nodenr: 3 coo: blade1 blade 1 root
+ 30 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:blade1 nodenr: 3 coo: blade1 blade 1 root
+ 31 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:blade1 nodenr: 10 coo: local blade 1 50% local e coo
+ 32 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:blade1 nodenr: 10 coo: local blade 1 50% local e coo
+ 33 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:blade1 nodenr: 10 coo: local blade 1 50% local e coo
+ 34 1.0 0.0 0.00 1.0 kN Fx coo: Force Fx Mbdy:blade1 nodenr: 1 coo: blade1 blade 1 root
+ 35 1.0 0.0 0.00 1.0 kN Fy coo: Force Fy Mbdy:blade1 nodenr: 1 coo: blade1 blade 1 root
+ 36 1.0 0.0 0.00 1.0 kN Fz coo: Force Fz Mbdy:blade1 nodenr: 1 coo: blade1 blade 1 root
+ 37 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:blade2 nodenr: 1 coo: blade2 blade 2 root
+ 38 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:blade2 nodenr: 1 coo: blade2 blade 2 root
+ 39 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:blade2 nodenr: 1 coo: blade2 blade 2 root
+ 40 1.0 0.0 0.00 1.0 kNm Mx coo: MomentMx Mbdy:blade3 nodenr: 1 coo: blade3 blade 3 root
+ 41 1.0 0.0 0.00 1.0 kNm My coo: MomentMy Mbdy:blade3 nodenr: 1 coo: blade3 blade 3 root
+ 42 1.0 0.0 0.00 1.0 kNm Mz coo: MomentMz Mbdy:blade3 nodenr: 1 coo: blade3 blade 3 root
+ 43 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:towertop E-nr: 1 Z-rel:1.00 coo: global tower top flange position
+ 44 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:towertop E-nr: 1 Z-rel:1.00 coo: global tower top flange position
+ 45 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:towertop E-nr: 1 Z-rel:1.00 coo: global tower top flange position
+ 46 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: blade1 blade 1 tip pos
+ 47 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: blade1 blade 1 tip pos
+ 48 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: blade1 blade 1 tip pos
+ 49 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: blade2 blade 2 tip pos
+ 50 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: blade2 blade 2 tip pos
+ 51 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: blade2 blade 2 tip pos
+ 52 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: blade3 blade 3 tip pos
+ 53 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: blade3 blade 3 tip pos
+ 54 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: blade3 blade 3 tip pos
+ 55 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: global blade 1 tip pos
+ 56 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: global blade 1 tip pos
+ 57 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: global blade 1 tip pos
+ 58 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: hub1 blade 1 tip pos hub coo
+ 59 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: hub1 blade 1 tip pos hub coo
+ 60 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade1 E-nr: 18 Z-rel:1.00 coo: hub1 blade 1 tip pos hub coo
+ 61 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: hub2 blade 2 tip pos hub coo
+ 62 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: hub2 blade 2 tip pos hub coo
+ 63 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade2 E-nr: 18 Z-rel:1.00 coo: hub2 blade 2 tip pos hub coo
+ 64 1.0 0.0 0.00 1.0 m State p State pos x Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: hub3 blade 3 tip pos hub coo
+ 65 1.0 0.0 0.00 1.0 m State p State pos y Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: hub3 blade 3 tip pos hub coo
+ 66 1.0 0.0 0.00 1.0 m State p State pos z Mbdy:blade3 E-nr: 18 Z-rel:1.00 coo: hub3 blade 3 tip pos hub coo
+ 67 1.0 0.0 0.00 1.0 rad/s omega t State_rot omega tz Mbdy:shaft E-nr: 4 Z-rel:1.00 coo: shaft
diff --git a/wetb/flex/tests/test_flex_io.py b/wetb/flex/tests/test_flex_io.py
new file mode 100644
index 0000000000000000000000000000000000000000..b82a693c0a8297d6aabd54d5193b38c7f2a9cb02
--- /dev/null
+++ b/wetb/flex/tests/test_flex_io.py
@@ -0,0 +1,27 @@
+'''
+Created on 11. apr. 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+from wetb import flex
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
+class Test(unittest.TestCase):
+
+
+ def test_load(self):
+ time, data, info = flex.load(tfp+"test1/test.int")
+ self.assertEqual(data.shape, (800,7))
+ self.assertEqual(info['attribute_names'][0], "WSP_gl._")
+ self.assertEqual(info['attribute_units'][0], "m/s")
+ self.assertEqual(info['attribute_descriptions'][0], "Free wind speed Vy, gl. coo, of gl. pos 0.75, 0.00, -40.75")
+
+ self.assertAlmostEqual(data[0, 1], 12.037,3)
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testName']
+ unittest.main()
\ No newline at end of file
diff --git a/wetb/flex/tests/test_sensor.py b/wetb/flex/tests/test_sensor.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a3e3ffeebe9f02e8dd34034dcb4469511afe235
--- /dev/null
+++ b/wetb/flex/tests/test_sensor.py
@@ -0,0 +1,33 @@
+'''
+Created on 9. sep. 2016
+
+@author: mmpe
+'''
+import os
+import unittest
+from wetb.flex import read_sensor_info
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
+class Test(unittest.TestCase):
+
+
+ def test_sensor_load(self):
+ sensor_info = read_sensor_info(tfp + "test_sensor_info/sensor")
+ nr, name, unit, description, _, _ = sensor_info[17]
+ self.assertEqual(nr, 18)
+ self.assertEqual(name, "Mz coo:")
+ self.assertEqual(unit, "kNm")
+ self.assertEqual(description, "MomentMz Mbdy:tower nodenr: 1 coo: tower tower base flange")
+
+ def test_sensor_load_name_stop(self):
+ sensor_info = read_sensor_info(tfp + "test_sensor_info/sensor"," ")
+ nr, name, unit, description, _, _ = sensor_info[17]
+ self.assertEqual(nr, 18)
+ self.assertEqual(name, "Mz")
+ self.assertEqual(unit, "kNm")
+ self.assertEqual(description, "coo: MomentMz Mbdy:tower nodenr: 1 coo: tower tower base flange")
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.test_sensor_load']
+ unittest.main()
diff --git a/wetb/gtsdf/__init__.py b/wetb/gtsdf/__init__.py
index 0c84199f52b04b0e90da6e7ac32a120f6e28a62b..4d48fa17a46720365d1de4717a393db9771bf8a9 100644
--- a/wetb/gtsdf/__init__.py
+++ b/wetb/gtsdf/__init__.py
@@ -67,7 +67,9 @@ class Dataset(object):
for i, n in enumerate(self.info['attribute_names']):
print (i,n)
raise e
-
+
+ def __contains__(self, name):
+ return name in self.info['attribute_names']
__all__ = sorted([m for m in set(dir()) - set(d)])
diff --git a/wetb/hawc2/ae_file.py b/wetb/hawc2/ae_file.py
index cee13d9f107fea4ad249ad514e480e70d61d2233..d3628fee19eebc8bc152588b7c572a703e01c179 100644
--- a/wetb/hawc2/ae_file.py
+++ b/wetb/hawc2/ae_file.py
@@ -45,14 +45,24 @@ class AEFile(object):
def _value(self, radius, column, set_nr=1):
ae_data = self.ae_sets[set_nr]
- return np.interp(radius, ae_data[:, 0], ae_data[:, column])
+ if radius:
+ return np.interp(radius, ae_data[:, 0], ae_data[:, column])
+ else:
+ return ae_data[:,column]
- def chord(self, radius, set_nr=1):
+ def chord(self, radius=None, set_nr=1):
return self._value(radius, 1, set_nr)
- def thickness(self, radius, set_nr=1):
+ def thickness(self, radius=None, set_nr=1):
return self._value(radius, 2, set_nr)
-
+
+ def radius_ae(self, radius=None, set_nr=1):
+ radii = self.ae_sets[set_nr][:,0]
+ if radius:
+ return radii[np.argmin(np.abs(radii-radius))]
+ else:
+ return radii
+
def pc_set_nr(self, radius, set_nr=1):
ae_data = self.ae_sets[set_nr]
index = np.searchsorted(ae_data[:, 0], radius)
@@ -68,6 +78,7 @@ class AEFile(object):
if __name__ == "__main__":
ae = AEFile(r"tests/test_files/NREL_5MW_ae.txt")
- print (ae.thickness(36))
+ print (ae.radius_ae(36))
+ print (ae.thickness())
print (ae.chord(36))
print (ae.pc_set_nr(36))
diff --git a/wetb/hawc2/at_time_file.py b/wetb/hawc2/at_time_file.py
index 10ed8583dd11b1bb280db37a1aa216d43bfcbde1..c84c75e6fdd38bb6646aac366556b261e40744fa 100644
--- a/wetb/hawc2/at_time_file.py
+++ b/wetb/hawc2/at_time_file.py
@@ -16,22 +16,37 @@ import numpy as np
class AtTimeFile(object):
"""Loads an at time file generated by HAWC2
- Note that the radius in this context is the curved radius
-
- >>> atfile = AtTimeFile("at_time.dat") # load file
- >>> atfile.attribute_names # Attribute names
+ Functions are generated for each sensor in the file from the template:
+ xxx(radius=None, curved_length=None)
+
+ E.g.
+ twist(radius=36) # Twist at radius 36m (Straight line radius, blade_radius must be specified in constructor)
+ twist(curved_length=36) # Twist 36m from the root along the (curved) center line
+ twist() # Twist at all aerodynamic calculation points
+
+
+ >>> at = AtTimeFile("at.dat", 86.3655) # load file
+ >>> at.attribute_names # Attribute names
['radius_s', 'twist', 'chord']
- >>> atfile[:3,1]) # first 3 twist rows
- [ 0. -0.775186 -2.91652 ]
- >>> atfile.twist()[:3]) # first 3 twist rows
- [ 0. -0.775186 -2.91652 ]
- >>> atfile.twist(curved_length=10) # Twist at curved_length = 10 (interpolated)
- -5.34743208242399
- >>> atfile.twist(radius=10) # Twist at curved_length = 10 (interpolated)
- -5.34743208242399
+ >>> at[:3,1]) # first 3 twist rows
+ [-14.5 -14.5002 -14.5007]
+ >>> at.twist()[:3]) # Twist first 3 twist rows
+ [-14.5 -14.5002 -14.5007]
+ >>> at.twist(curved_length=36) # Twist at curved_length = 10 (interpolated)
+ -5.172195702789108
+ >>> at.twist(radius=36) # Twist at curved_length = 10 (interpolated)
+ -5.162084567646019
"""
- def __init__(self, filename, blade_radius=None):
- self.blade_radius = blade_radius
+ def __init__(self, filename, bladetip_radius=None):
+ """
+ Parameters
+ ----------
+ filename : string
+ Filename
+ bladetip_radius : int, float or None, optional
+ Radius of blade tip. Used to convert from curved length to radius
+ """
+ self.blade_radius = bladetip_radius
with open(filename, encoding='utf-8') as fid:
lines = fid.readlines()
self.attribute_names = lines[2].lower().replace("#", "").split()
@@ -43,7 +58,7 @@ class AtTimeFile(object):
if radius is None and curved_length is None:
return self.data[:, column]
elif radius is not None:
- assert self.blade_radius is not None, "blade_radius must be specified in __init__ when requesting value of radius (alternatively you can request for curved_length)"
+ assert self.blade_radius is not None, "bladetip_radius must be specified in __init__ when requesting value of radius (alternatively you can request for curved_length)"
return np.interp(radius/self.blade_radius, self.data[:, 0]/self.data[-1, 0], self.data[:, column])
else:
return np.interp(curved_length, self.data[:, 0], self.data[:, column])
@@ -51,7 +66,7 @@ class AtTimeFile(object):
for column, att_name in enumerate(self.attribute_names):
setattr(self, att_name, func_factory(column))
- def ac_radius(self, radius=None):
+ def radius_ac(self, radius=None):
"""Radius (curved distance) of aerodynamic calculation point(s)
Parameters
@@ -66,7 +81,7 @@ class AtTimeFile(object):
Radius of calculation points or radius of calculation point nearest radius
"""
if radius is None:
- return self.radius_s(radius)
+ return self.radius_s()
else:
return self.radius_s()[np.argmin(np.abs(self.radius_s() - radius))]
@@ -78,8 +93,18 @@ class AtTimeFile(object):
if __name__ == "__main__":
- at = AtTimeFile(r"tests/test_files/at_time.dat")
- print (at.attribute_names)
- print (at.twist(36))
- print (at.chord(36))
+ at = AtTimeFile(r"tests/test_files/at.dat", 86.3655) # load file
+ print (at.attribute_names) # Attribute names
+ print (at[:3,1]) # first 3 twist rows
+ print (at.twist()[:3]) # first 3 twist rows
+ print (at.twist(curved_length=36)) # Twist at curved_length = 10 (interpolated)
+ print (at.twist(radius=36)) # Twist at curved_length = 10 (interpolated)
+
+# at = AtTimeFile(r"tests/test_files/at.dat", 86.3655)
+# print (at.attribute_names)
+# print (at.radius_s())
+# print (at.twist(curved_length=36))
+# print (at.chord(curved_length=36))
+# print (at.twist(curved_length=36))
+# print (at.twist(radius=36))
diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/bladeinfo.py
new file mode 100644
index 0000000000000000000000000000000000000000..745685f31f1550761127592954c20890b2f658dc
--- /dev/null
+++ b/wetb/hawc2/bladeinfo.py
@@ -0,0 +1,210 @@
+'''
+Created on 3. maj 2017
+
+@author: mmpe
+'''
+import numpy as np
+from wetb.hawc2.at_time_file import AtTimeFile
+from wetb.hawc2.pc_file import PCFile
+from wetb.hawc2.htc_file import HTCFile
+from scipy.interpolate.interpolate import interp1d
+import os
+from wetb.utils.geometry import rad
+from wetb.hawc2.st_file import StFile
+
+
+
+class BladeInfo(object):
+ def twist(self, radius=None):
+ """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
+
+ Parameters
+ ---------
+ radius : float or array_like, optional
+ radius [m] of interest
+ If None (default) the twist of all points are returned
+
+ Returns
+ -------
+ twist [deg] : float
+ """
+ raise NotImplementedError()
+
+ def chord(self, radius=None):
+ """Aerodynamic chord
+
+ Parameters
+ ---------
+ radius : float or array_like, optional
+ radius [m] of interest
+ If None (default) the twist of all points are returned
+
+ Returns
+ -------
+ chord [m] : float
+ """
+ raise NotImplementedError()
+
+
+
+ def c2nd(self, radius_nd):
+ """Return center line position
+
+ Parameters
+ ---------
+ radius_nd : float or array_like, optional
+ non dimensional radius
+
+ Returns
+ -------
+ x,y,z : float
+ """
+
+class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
+ """Provide HAWC2 info about a blade
+
+ From AE file:
+ - chord(radius=None, set_nr=1):
+ - thickness(radius=None, set_nr=1)
+ - radius_ae(radius=None, set_nr=1)
+
+ From PC file
+ - CL(radius, alpha, ae_set_nr=1)
+ - CD(radius, alpha, ae_set_nr=1)
+ - CM(radius, alpha, ae_set_nr=1)
+
+ From at_time_filename
+ - attribute_names
+ - xxx(radius=None, curved_length=None) # xxx for each attribute name
+ - radius_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
+
+ From ST file
+ - radius_st(radius=None, mset=1, set=1)
+ - xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
+ with template
+
+ """
+
+ def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+
+ if isinstance(htcfile, str):
+ assert htcfile.lower().endswith('.htc')
+ htcfile = HTCFile(htcfile)
+
+ blade_name = blade_name or htcfile.aero.link[2]
+ s = htcfile.new_htc_structure
+ at_time_filename = at_time_filename or os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat")
+ pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
+ ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
+
+ mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+ mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
+ st_filename = st_filename or os.path.join(htcfile.modelpath, mainbody_blade.timoschenko_input.filename[0])
+
+ if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
+ PCFile.__init__(self, pc_filename, ae_filename)
+ blade_radius = self.ae_sets[1][-1,0]
+
+ if os.path.isfile(st_filename):
+ StFile.__init__(self, st_filename)
+ if os.path.isfile(at_time_filename):
+ AtTimeFile.__init__(self, at_time_filename, blade_radius)
+
+
+ self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
+
+ #self.c2nd = lambda x : interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')(np.max([np.min([x, np.ones_like(x)], 0), np.zeros_like(x) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+ x, y, z, twist = self.hawc2_splines()
+ def interpolate(r):
+ r = (np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+ return np.array([np.interp(r, np.array(z) / z[-1], xyz) for xyz in [x,y,z, twist]]).T
+
+ self.c2nd = interpolate #lambda r : interp1d(np.array(z) / z[-1], np.array([x, y, z]).T, axis=0, kind=1)(np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+
+ def hawc2_splines(self):
+ curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+ curve_z_nd = curve_z / curve_z[-1]
+
+ def akima(x, y):
+ n = len(x)
+ var = np.zeros((n + 3))
+ z = np.zeros((n))
+ co = np.zeros((n, 4))
+ for i in range(n - 1):
+ var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
+ var[n + 1] = 2 * var[n] - var[n - 1]
+ var[n + 2] = 2 * var[n + 1] - var[n]
+ var[1] = 2 * var[2] - var[3]
+ var[0] = 2 * var[1] - var[2]
+
+ for i in range(n):
+ wi1 = abs(var[i + 3] - var[i + 2])
+ wi = abs(var[i + 1] - var[i])
+ if (wi1 + wi) == 0:
+ z[i] = (var[i + 2] + var[i + 1]) / 2
+ else:
+ z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
+
+ for i in range(n - 1):
+ dx = x[i + 1] - x[i]
+ a = (z[i + 1] - z[i]) * dx
+ b = y[i + 1] - y[i] - z[i] * dx
+ co[i, 0] = y[i]
+ co[i, 1] = z[i]
+ co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
+ co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
+ co[n - 1, 0] = y[n - 1]
+ co[n - 1, 1] = z[n - 1]
+ co[n - 1, 2] = 0
+ co[n - 1, 3] = 0
+ return co
+
+ def coef2spline(s, co):
+ x, y = [], []
+ for i, c in enumerate(co.tolist()[:-1]):
+ p = np.poly1d(c[::-1])
+ z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2)
+ x.extend(s[i] + z)
+ y.extend(p(z))
+ return y
+
+ x, y, z, twist = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(4)]
+ return x, y, z, twist
+
+ def c2def_twist(self, radius=None):
+ if radius is None:
+ return self.c2def[:, 3]
+ else:
+ return np.interp(radius, self.c2def[:, 2], self.c2def[:, 3])
+
+
+
+class H2aeroBladeInfo(H2BladeInfo):
+
+ def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
+ """
+ at_time_filename: file name of at time file containing twist and chord data
+ """
+ PCFile.__init__(self, pc_filename, ae_filename)
+ blade_radius = self.ae_sets[1][-1,0]
+ AtTimeFile.__init__(self, at_time_filename, blade_radius)
+
+ assert('twist' in self.attribute_names)
+ htcfile = HTCFile(htc_filename)
+
+
+ self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
+ #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
+
+ ac_radii = self.ac_radius()
+ c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
+ self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
+ #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
+
+ def c2nd(self, r_nd):
+ r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
+ r_nd_max = np.ones_like(r_nd)
+ r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
+ return self._c2nd(r_nd)
+
+
diff --git a/wetb/hawc2/htc_contents.py b/wetb/hawc2/htc_contents.py
index 5337613d47c902590c483da4ec6fd8e22d501596..93b58f077cfd2fb2cdb12c9cfdb663f0a813cb76 100644
--- a/wetb/hawc2/htc_contents.py
+++ b/wetb/hawc2/htc_contents.py
@@ -120,7 +120,7 @@ class HTCContents(object):
self._add_contents(section)
return section
- def add_line(self, name, values, comments):
+ def add_line(self, name, values, comments=""):
line = HTCLine(name, values, comments)
self._add_contents(line)
return line
diff --git a/wetb/hawc2/pc_file.py b/wetb/hawc2/pc_file.py
index f42cc1076518a42f14b38cb307b65c6847a28a39..bd0da0a7403545b893cbea7ac1a2b4dd5f432689 100644
--- a/wetb/hawc2/pc_file.py
+++ b/wetb/hawc2/pc_file.py
@@ -43,7 +43,7 @@ class PCFile(AEFile):
with open (filename) as fid:
lines = fid.readlines()
nsets = int(lines[0].split()[0])
- self.sets = {}
+ self.pc_sets = {}
lptr = 1
for nset in range(1, nsets + 1):
nprofiles = int(lines[lptr].split()[0])
@@ -59,12 +59,12 @@ class PCFile(AEFile):
thicknesses.append(thickness)
profiles.append(data)
lptr += n_rows
- self.sets[nset] = (np.array(thicknesses), profiles)
+ self.pc_sets[nset] = (np.array(thicknesses), profiles)
def _Cxxx(self, radius, alpha, column, ae_set_nr=1):
thickness = self.thickness(radius, ae_set_nr)
pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
- thicknesses, profiles = self.sets[pc_set_nr]
+ thicknesses, profiles = self.pc_sets[pc_set_nr]
index = np.searchsorted(thicknesses, thickness)
if index == 0:
index = 1
@@ -112,7 +112,6 @@ class PCFile(AEFile):
return self._Cxxx(radius, alpha, 2, ae_set_nr)
def CM(self, radius, alpha, ae_set_nr=1):
-
return self._Cxxx(radius, alpha, 3, ae_set_nr)
if __name__ == "__main__":
diff --git a/wetb/hawc2/st_file.py b/wetb/hawc2/st_file.py
index f03efafb96d5b93b58c330e427abaed48f09e2cb..7859144fab6a442ff110c681297e46fd4a32690a 100644
--- a/wetb/hawc2/st_file.py
+++ b/wetb/hawc2/st_file.py
@@ -68,8 +68,9 @@ class StFile(object):
def __init__(self, filename):
with open (filename) as fid:
txt = fid.read()
- no_maindata_sets = int(txt.replace("#","").strip()[0])
- assert no_maindata_sets == txt.count("#")
+# Some files starts with first set ("#1...") with out specifying number of sets
+# no_maindata_sets = int(txt.strip()[0])
+# assert no_maindata_sets == txt.count("#")
self.main_data_sets = {}
for mset in txt.split("#")[1:]:
mset_nr = int(mset.strip().split()[0])
@@ -91,7 +92,7 @@ class StFile(object):
radius = self.radius(None, mset_nr, set_nr)
return np.interp(radius, st_data[:, 0], st_data[:, column])
- def radius(self, radius=None, mset=1, set=1):
+ def radius_st(self, radius=None, mset=1, set=1):
r = self.main_data_sets[mset][set][:, 0]
if radius is None:
return r
diff --git a/wetb/hawc2/tests/test_AtTimeFile.py b/wetb/hawc2/tests/test_AtTimeFile.py
index 6f17562bab8741498c04106893e0541d9f50a17f..09573fd33ba82feb1f5dc5d52cc5ba1ec1ea0b21 100644
--- a/wetb/hawc2/tests/test_AtTimeFile.py
+++ b/wetb/hawc2/tests/test_AtTimeFile.py
@@ -15,7 +15,7 @@ import numpy as np
import os
-
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
class TestAtTimeFile(unittest.TestCase):
def setUp(self):
@@ -24,7 +24,7 @@ class TestAtTimeFile(unittest.TestCase):
def test_doc_examples(self):
- atfile = AtTimeFile(self.testfilepath + "at_time.dat", blade_radius=20.501) # load file
+ atfile = AtTimeFile(self.testfilepath + "at_time.dat", bladetip_radius=20.501) # load file
self.assertEqual(atfile.attribute_names, ['radius_s', 'twist', 'chord'])
np.testing.assert_array_equal(atfile[:3, 1], [ 0., -0.775186, -2.91652 ])
np.testing.assert_array_equal(atfile.twist()[:3], [ 0. , -0.775186 , -2.91652 ])
@@ -48,7 +48,7 @@ class TestAtTimeFile(unittest.TestCase):
self.assertEqual(atfile.chord(curved_length=9), 1.3888996578373045)
def test_at_time_file_at_radius(self):
- atfile = AtTimeFile(self.testfilepath + "at_time.dat", blade_radius=20.501/2)
+ atfile = AtTimeFile(self.testfilepath + "at_time.dat", bladetip_radius=20.501/2)
self.assertEqual(atfile.radius_s(radius=9/2), 9)
self.assertEqual(atfile.twist(radius=9/2), -6.635983309665461)
self.assertEqual(atfile.chord(radius=9/2), 1.3888996578373045)
@@ -56,9 +56,9 @@ class TestAtTimeFile(unittest.TestCase):
def test_at_time_file_radius(self):
atfile = AtTimeFile(self.testfilepath + "at_time.dat")
- self.assertEqual(atfile.ac_radius()[12], 10.2505)
- self.assertEqual(atfile.ac_radius(10), 10.2505)
- self.assertEqual(atfile.ac_radius(10.5), 10.2505)
+ self.assertEqual(atfile.radius_ac()[12], 10.2505)
+ self.assertEqual(atfile.radius_ac(10), 10.2505)
+ self.assertEqual(atfile.radius_ac(10.5), 10.2505)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
diff --git a/wetb/hawc2/tests/test_files/at.dat b/wetb/hawc2/tests/test_files/at.dat
new file mode 100644
index 0000000000000000000000000000000000000000..f93af10f4f3c18596f1dd613503a5bc82847fffe
--- /dev/null
+++ b/wetb/hawc2/tests/test_files/at.dat
@@ -0,0 +1,53 @@
+ # Version ID : HAWC2MB 12.4
+ # Radial output at time: 0.120000000000000
+ # Radius_s # twist # Chord
+ 0.00000E+00 -1.45000E+01 5.38000E+00
+ 8.88594E-02 -1.45002E+01 5.38000E+00
+ 3.55072E-01 -1.45007E+01 5.38000E+00
+ 7.97544E-01 -1.45012E+01 5.38000E+00
+ 1.41446E+00 -1.45013E+01 5.38000E+00
+ 2.20328E+00 -1.45008E+01 5.38000E+00
+ 3.16076E+00 -1.44998E+01 5.38000E+00
+ 4.28297E+00 -1.44959E+01 5.38000E+00
+ 5.56530E+00 -1.44880E+01 5.38150E+00
+ 7.00248E+00 -1.44609E+01 5.40757E+00
+ 8.58860E+00 -1.43497E+01 5.47316E+00
+ 1.03171E+01 -1.40434E+01 5.57870E+00
+ 1.21810E+01 -1.34021E+01 5.71677E+00
+ 1.41726E+01 -1.24242E+01 5.87200E+00
+ 1.62836E+01 -1.11781E+01 6.02007E+00
+ 1.85054E+01 -9.82623E+00 6.13078E+00
+ 2.08289E+01 -8.68972E+00 6.18999E+00
+ 2.32446E+01 -7.84255E+00 6.19714E+00
+ 2.57424E+01 -7.19893E+00 6.15543E+00
+ 2.83122E+01 -6.63223E+00 6.06940E+00
+ 3.09433E+01 -6.10664E+00 5.94438E+00
+ 3.36250E+01 -5.60732E+00 5.78666E+00
+ 3.63463E+01 -5.10875E+00 5.60337E+00
+ 3.90959E+01 -4.60713E+00 5.40154E+00
+ 4.18626E+01 -4.07749E+00 5.18393E+00
+ 4.46350E+01 -3.51118E+00 4.95580E+00
+ 4.74017E+01 -2.93766E+00 4.72176E+00
+ 5.01513E+01 -2.36236E+00 4.48569E+00
+ 5.28726E+01 -1.79469E+00 4.25127E+00
+ 5.55543E+01 -1.25288E+00 4.02160E+00
+ 5.81854E+01 -7.42980E-01 3.79894E+00
+ 6.07552E+01 -2.63044E-01 3.58510E+00
+ 6.32530E+01 1.76184E-01 3.38191E+00
+ 6.56687E+01 5.69023E-01 3.19012E+00
+ 6.79922E+01 9.24087E-01 3.01062E+00
+ 7.02140E+01 1.25492E+00 2.84386E+00
+ 7.23251E+01 1.55945E+00 2.68993E+00
+ 7.43167E+01 1.84168E+00 2.54872E+00
+ 7.61805E+01 2.10535E+00 2.42040E+00
+ 7.79091E+01 2.34909E+00 2.29766E+00
+ 7.94952E+01 2.57177E+00 2.17211E+00
+ 8.09324E+01 2.77156E+00 2.04065E+00
+ 8.22147E+01 2.94258E+00 1.90297E+00
+ 8.33369E+01 3.08259E+00 1.75760E+00
+ 8.42945E+01 3.19421E+00 1.59874E+00
+ 8.50832E+01 3.28140E+00 1.41325E+00
+ 8.56999E+01 3.34684E+00 1.20310E+00
+ 8.61425E+01 3.39238E+00 1.00803E+00
+ 8.64090E+01 3.41914E+00 8.77215E-01
+ 8.64979E+01 3.42796E+00 8.33540E-01
diff --git a/wetb/hawc2/tests/test_st_file.py b/wetb/hawc2/tests/test_st_file.py
index 8d794c6ea2d6258f23d713ffe94c5d1be3856592..8b678b65b1a466d61749dd380938d515a0c18ad8 100644
--- a/wetb/hawc2/tests/test_st_file.py
+++ b/wetb/hawc2/tests/test_st_file.py
@@ -19,8 +19,8 @@ class TestStFile(unittest.TestCase):
def test_stfile(self):
st = StFile(testfilepath + "DTU_10MW_RWT_Blade_st.dat")
- self.assertEqual(st.radius()[2], 3.74238)
- self.assertEqual(st.radius(3), 3.74238)
+ self.assertEqual(st.radius_st()[2], 3.74238)
+ self.assertEqual(st.radius_st(3), 3.74238)
self.assertEqual(st.x_e(67.7351), 4.4320990737400E-01)
self.assertEqual(st.E(3.74238, 1, 1), 1.2511695058500E+10)
self.assertEqual(st.E(3.74238, 1, 2), 1.2511695058500E+27)
diff --git a/wetb/signal/fix/_rotor_position.py b/wetb/signal/fix/_rotor_position.py
index 22aac76de54ba82a0cd4d21f501885e9a7364cb9..3876204dc3c039b8bfa047b028f489f0e8053107 100644
--- a/wetb/signal/fix/_rotor_position.py
+++ b/wetb/signal/fix/_rotor_position.py
@@ -5,6 +5,8 @@ Created on 30. mar. 2017
'''
import numpy as np
from wetb.signal.fit._spline_fit import spline_fit
+from wetb.signal.filters._differentiation import differentiation
+
def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
"""Rotor position fitted with spline
@@ -33,7 +35,7 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
from wetb.signal.subset_mean import revolution_trigger
t = np.arange(len(rotor_position))
- indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_no_round_diff=4)
+ indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=4)
rp = rotor_position[:].copy()
@@ -69,7 +71,9 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
plt.plot(t/sample_frq, spline_fit(x,y)(t)-t*a, label='Spline (detrended)')
plt.legend()
plt.show()
- return spline_fit(x,y)(t)%360
+ fit = spline_fit(x,y)(t)
+ fit[differentiation(fit, "left")<0]=np.nan
+ return fit%360
def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
"""Find the optimal fix_dt parameter for fix_rotor_position (function above).
@@ -95,7 +99,8 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
"""
from wetb.signal.filters import differentiation
def err(i):
- rpm_pos = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))%180 / 360 * sample_frq * 60
+ drp = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))
+ rpm_pos = drp%180 / 360 * sample_frq * 60
return np.sum((rpm_pos - rotor_speed)**2)
best = 27
@@ -110,4 +115,69 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
plt.show()
return best
-
\ No newline at end of file
+
+def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
+ """Rotor position fitted with spline
+
+ Parameters
+ ----------
+ rotor_position : array_like
+ Rotor position [deg] (0-360)
+ sample_frq : int or float
+ Sample frequency [Hz]
+ rotor_speed : array_like
+ Rotor speed [RPM]
+ fix_dt : int, float or None, optional
+ Time distance [s] between spline fix points\n
+ If None (default) a range of seconds is tested and the result that minimize the RMS
+ between differentiated rotor position fit and rotor speed is used.\n
+ Note that a significant speed up is achievable by specifying the parameter
+ plt : PyPlot or None
+ If PyPlot a visual interpretation is plotted
+
+ Returns
+ -------
+ y : nd_array
+ Fitted rotor position
+ """
+
+ from wetb.signal.subset_mean import revolution_trigger
+
+ t = np.arange(len(rotor_position))/sample_frq
+ indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=max_rev_diff)
+
+ rotor_position = rotor_position[:]%360
+
+ rp_fit = fix_rotor_position(rotor_position, sample_frq, rotor_speed, None)
+
+ if plt:
+ #a = (rp.max()-rp.min())/t.max()
+ #plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
+ f, axarr = plt.subplots(2)
+ print (rp_fit)
+ axarr[0].plot(t, rotor_position)
+ axarr[0].plot(indexes/sample_frq, rotor_position[indexes],'.')
+ axarr[0].plot(t, rp_fit)
+ axarr[1].plot(t, rotor_speed)
+ axarr[1].plot(indexes/sample_frq,60/( differentiation(indexes)/sample_frq))
+ print (t[170:200])
+ drp = differentiation(rp_fit)
+ #drp[(drp<0)]= 0
+ axarr[1].plot(t, drp%180 / 360 * sample_frq * 60, label='fix')
+ plt.legend()
+
+ i1,i2 = indexes[0], indexes[-1]
+ print ("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2])%180)/360)
+ print ("Rev from rotor speed", np.mean(rotor_speed[i1:i2])*(i2-i1)/60/sample_frq)
+ print ("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360)
+
+ print ("Mean RPM from rotor speed", np.mean(rotor_speed))
+ print ("Mean RPM from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360 / ((i2-i1)/60/sample_frq))
+
+ spr1 = (np.diff(indexes)/sample_frq)
+
+ #rs1 = 60/( np.diff(indexes)/sample_frq)
+ spr2 = np.array([60/rotor_speed[i1:i2].mean() for i1,i2 in zip(indexes[:-1], indexes[1:])])
+
+ err = spr1-spr2
+ print (err.max())
\ No newline at end of file
diff --git a/wetb/signal/subset_mean.py b/wetb/signal/subset_mean.py
index 0de3ea896f8185e932998ddd85229b97523e2a9a..7578c11ac477856a237b7ec00ea0f2e418ba6695 100644
--- a/wetb/signal/subset_mean.py
+++ b/wetb/signal/subset_mean.py
@@ -9,6 +9,7 @@ import unittest
from wetb.utils.geometry import rpm2rads
from _collections import deque
from tables.tests.test_index_backcompat import Indexes2_0TestCase
+from wetb.signal.filters._differentiation import differentiation
def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
@@ -136,7 +137,7 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance=
else:
return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step]
-def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_diff=1):
+def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1):
"""Returns one index per revolution (minimum rotor position)
Parameters
@@ -155,16 +156,12 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
if isinstance(rotor_speed, (float, int)):
rotor_speed = np.ones_like(rotor_position)*rotor_speed
deg_per_sample = rotor_speed*360/60/sample_frq
- sample_per_round = 1/(rotor_speed/60/sample_frq)
- thresshold = deg_per_sample.max()*2
-
- nround_rotor_speed = np.nansum(rotor_speed/60/sample_frq)
-
- mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+ thresshold = deg_per_sample.max()*3
+ drp = np.diff(rotor_position)%360
+ rp = rotor_position
+
- nround_rotor_position = np.nansum(np.diff(rotor_position)%mod)/360
- #assert abs(nround_rotor_position-nround_rotor_speed)<max_no_round_diff, "No of rounds from rotor_position (%.2f) mismatch with no_rounds from rotor_speed (%.2f)"%(nround_rotor_position, nround_rotor_speed)
- #print (nround_rotor_position, nround_rotor_speed)
+
rp = np.array(rotor_position).copy()
#filter degree increase > thresshold
@@ -176,10 +173,47 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
# Best lower is the first lower after upper
best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
upper2lower = best_lower - upper_indexes
- best_lower = best_lower[upper2lower<upper2lower.mean()*2]
- #max_dist_error = max([np.abs((i2-i1)- np.mean(sample_per_round[i1:i2])) for i1,i2 in zip(best_lower[:-1], best_lower[1:])])
- #assert max_dist_error < sample_frq/5, max_dist_error
- return best_lower
+ trigger_indexes = best_lower[upper2lower<upper2lower.mean()*2].tolist()
+
+ if len(trigger_indexes)>1:
+
+
+
+ rpm_rs = np.array([rev.mean() for rev in np.split(rotor_speed, trigger_indexes)[1:-1]])
+ rpm_i = 1/np.diff(trigger_indexes)*60*sample_frq
+ spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+ spr_i = np.diff(trigger_indexes)
+
+
+ while np.any(spr_rs>spr_i*1.9):
+ i = np.where(spr_rs>spr_i*1.9)[0][0]
+ if np.abs(spr_i[i-1] + spr_i[i] - spr_rs[i])< np.abs(spr_i[i] + spr_i[i+1] - spr_rs[i]):
+ del trigger_indexes[i]
+ else:
+ del trigger_indexes[i+1]
+ spr_i = np.diff(trigger_indexes)
+ spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+
+
+
+ # if a revolution is too long add trigger
+ if np.any(rpm_rs>rpm_i*2.1):
+ # >1 missing triggers
+ raise NotImplementedError
+ trigger_indexes.extend([np.mean(trigger_indexes[i:i+2]) for i in np.where(rpm_rs>rpm_i*1.9)[0]])
+ trigger_indexes = np.sort(trigger_indexes).astype(np.int)
+
+
+ i1,i2 = trigger_indexes[0], trigger_indexes[-1]
+ nround_rotor_speed = np.nansum(rotor_speed[i1:i2]/60/sample_frq)
+ #mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+ nround_rotor_position = len(trigger_indexes)-1 #np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
+ if max_rev_diff is not None:
+ diff_pct = abs(nround_rotor_position-nround_rotor_speed)/nround_rotor_position*100
+ assert diff_pct<max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%"%(nround_rotor_position, nround_rotor_speed, diff_pct)
+
+
+ return trigger_indexes
def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
diff --git a/wetb/signal/tests/test_fix.py b/wetb/signal/tests/test_fix.py
index 8663986285a6fdb9f3a37f487b751b619bd9cc83..86c9ce472931b9aacd3d5b0d081d0639988f6a2c 100644
--- a/wetb/signal/tests/test_fix.py
+++ b/wetb/signal/tests/test_fix.py
@@ -42,6 +42,8 @@ class TestFix(unittest.TestCase):
self.assertEqual(find_fix_dt(ds.azi, sample_frq, ds.Rot_cor), 4)
+
+
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testRotorPositionFix']
diff --git a/wetb/signal/tests/test_subset_mean.py b/wetb/signal/tests/test_subset_mean.py
index 23642ae7509bc80188499122807d5b22ddfce2fc..8c4ecd862d112fad16d5db3f85c86604d5e45ac8 100644
--- a/wetb/signal/tests/test_subset_mean.py
+++ b/wetb/signal/tests/test_subset_mean.py
@@ -92,21 +92,32 @@ class TestSubsetMean(unittest.TestCase):
x2 = np.random.randint(0, len(rotor_position),10)
print (list(x2))
else:
- x1 = [447, 854, 595, 804, 847, 488, 412, 199, 675, 766]
+ x1 = [447, 854, 847, 488, 412, 199, 675]
x2 = [92, 647, 821, 422, 33, 159, 369, 99, 157, 464]
rotor_position[x1] += 360
rotor_position[x2] -= 360
rotor_position[90] = 180
+ rotor_position[270:276] = 15
+ rotor_position[598:602] = [360,360,0,0]
+ rotor_position[748:752] = [360,360,0,0]
indexes = revolution_trigger(rotor_position, 20,1/(90/20/60))
- np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
if 0:
import matplotlib.pyplot as plt
plt.plot(rotor_position)
plt.plot(indexes, np.zeros_like(indexes),'.')
plt.show()
+ np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
+
+ def test_revolution_trigger_max_rev_diff(self):
+ rotor_position = np.arange(0.,360*10,4)
+ rotor_position += np.random.random(len(rotor_position))
+ rotor_position = rotor_position % 360
+ rotor_speed = 1/(90/20/60)
+ revolution_trigger(rotor_position, 20, rotor_speed*1.02, max_rev_diff=3)
+ self.assertRaisesRegex(AssertionError, "No of revolution mismatch", revolution_trigger, rotor_position, 20, rotor_speed*1.02, max_rev_diff=1)
+
-
if __name__ == "__main__":
unittest.main()
diff --git a/wetb/wind/__init__.py b/wetb/wind/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..ae17e9a567028038ec2029f25b2abf6d80ffa4ee 100644
--- a/wetb/wind/__init__.py
+++ b/wetb/wind/__init__.py
@@ -0,0 +1 @@
+from .turbulence import *
\ No newline at end of file
diff --git a/wetb/wind/turbulence/mann_parameters.py b/wetb/wind/turbulence/mann_parameters.py
index b6a1d806cf5d083cdd0510d391285bedfd7f922c..b0a94bc7c9dfa87895e986769c1f45013c030642 100644
--- a/wetb/wind/turbulence/mann_parameters.py
+++ b/wetb/wind/turbulence/mann_parameters.py
@@ -254,20 +254,18 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
return ae
-def plot_spectra(ae, L, G, k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, show=True, plt=None):
+def plot_spectra(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
# if plt is None:
# import matplotlib.pyplot as plt
_plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
plot_mann_spectra(ae, L, G, "-", mean_u, plt)
- if show:
- plt.show()
-def _plot_spectra(k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+def _plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
if plt is None:
import matplotlib.pyplot as plt
bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
def plot(xx, label, color, plt):
- plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color)
+ plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color, label=label)
plot(buu, 'uu', 'r', plt)
if (bvv) is not None:
plot(bvv, 'vv', 'g', plt)
@@ -280,8 +278,9 @@ def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv
import matplotlib.pyplot as plt
mf = 10 ** (np.linspace(-4, 1, 1000))
muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
+ plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G))
if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
- if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
+ if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)