diff --git a/docs/developer-guide.md b/docs/developer-guide.md
index e27a107fc94e027e175d585e97870972c16bab0c..f65d449a04820674bf4162acd306eea0af82a5c9 100644
--- a/docs/developer-guide.md
+++ b/docs/developer-guide.md
@@ -146,6 +146,8 @@ use ```deactivate``` to deactivate the environment.
 - [parimeko](http://www.paramiko.org/)
 - [sshtunnel](https://github.com/pahaz/sshtunnel)
 - [pandoc](http://pandoc.org/) , [pypandoc](https://pypi.python.org/pypi/pypandoc):
+- [click](https://click.palletsprojects.com/en/7.x/)
+- [jinja2](https://jinja.palletsprojects.com/en/2.10.x/)
 convert markdown formatted readme file to rst for PyPi compatibility. See also
 issue #22. ```pandoc``` is available in Anaconda. When installing
 ```pypandoc``` via pip, you have to install ```pandoc``` via your package
@@ -157,7 +159,7 @@ Install the necessary Python dependencies using the conda package manager:
 
 ```
 >> conda install setuptools_scm future h5py pytables pytest pytest-cov nose sphinx blosc pbr paramiko
->> conda install scipy pandas matplotlib cython xlrd coverage xlwt openpyxl psutil pandoc twine pypandoc
+>> conda install scipy pandas matplotlib cython xlrd coverage xlwt openpyxl psutil pandoc twine pypandoc click jinja2
 >> conda install -c conda-forge sshtunnel --no-deps
 ```
 
@@ -290,4 +292,3 @@ In case of problems:
 will fail. This means commit hashes can not be part of the version number.
 Note that when your git working directory is not clean, the scheme for automatic
 versioning number will add ```dirty``` to the version number.
-
diff --git a/setup.py b/setup.py
index 9902b35baba1648a54b56e6e3138a08343af1984..08f404f0a677461c6dcc5cba3de8871d9c0367e6 100644
--- a/setup.py
+++ b/setup.py
@@ -53,7 +53,9 @@ def setup_package():
                         'openpyxl',
                         'psutil',
                         'six',
-                        'sshtunnel']
+                        'sshtunnel',
+                        'click',
+                        'jinja2',]
     build_requires = ['cython']
     setup(install_requires=install_requires,
           setup_requires=install_requires + build_requires + sphinx,
diff --git a/wetb/dlb/tests/test_iec64100_1_dlb.py b/wetb/dlb/tests/test_iec64100_1_dlb.py
index 0e8ed29d7f99c65cf87d340c1f14ebc6df8cdcd9..f151d204d2453fdbe6201e832d873e25d3b67966 100644
--- a/wetb/dlb/tests/test_iec64100_1_dlb.py
+++ b/wetb/dlb/tests/test_iec64100_1_dlb.py
@@ -45,6 +45,7 @@ def test_DLC12(writer):
     dlc12 = DTU_IEC64100_1_Ref_DLB(iec_wt_class='1A', Vin=4, Vout=26, Vr=10, D=180, z_hub=90)['DLC12']
     assert len(dlc12) == 216  # 12 wsp, 3 wdir, 6 seeds
     writer.from_pandas(dlc12[::24][:2])
+    writer.write_all(path + "htc/DLC12")
     npt.assert_array_equal(sorted(os.listdir(path + "htc/DLC12")),
                            ['DLC12_wsp04_wdir350_s1001.htc', 'DLC12_wsp06_wdir000_s1101.htc'])
     htc = HTCFile(path + "htc/DLC12/DLC12_wsp04_wdir350_s1001.htc")
@@ -58,6 +59,7 @@ def test_DLC21(writer):
     dlc = DTU_IEC64100_1_Ref_DLB(iec_wt_class='1A', Vin=4, Vout=26, Vr=10, D=180, z_hub=90)['DLC21']
     assert len(dlc) == 144  # 12 wsp, 3 wdir, 4 seeds
     writer.from_pandas(dlc[::16][:2])
+    writer.write_all(path + "htc/DLC21")
     npt.assert_array_equal(sorted(os.listdir(path + "htc/DLC21")),
                            ['DLC21_wsp04_wdir350_s1001.htc', 'DLC21_wsp06_wdir000_s1101.htc'])
     htc = HTCFile(path + "htc/DLC21/DLC21_wsp04_wdir350_s1001.htc")
@@ -72,6 +74,7 @@ def test_DLC22y(writer):
     dlc = DTU_IEC64100_1_Ref_DLB(iec_wt_class='1A', Vin=4, Vout=26, Vr=10, D=180, z_hub=90)['DLC22y']
     assert len(dlc) == 276  # 12 wsp, 23 wdir, 1 seeds
     writer.from_pandas(dlc[::24][:2])
+    writer.write_all(path + "htc/DLC22y")
     npt.assert_array_equal(sorted(os.listdir(path + "htc/DLC22y")),
                            ['DLC22y_wsp04_wdir015_s1001.htc', 'DLC22y_wsp06_wdir030_s1101.htc'])
     htc = HTCFile(path + "htc/DLC22y/DLC22y_wsp04_wdir015_s1001.htc")
diff --git a/wetb/hawc2/hawc2_input_writer.py b/wetb/hawc2/hawc2_input_writer.py
index 4accbd6a93c269b9b8358c44f69439b1851b7b9a..3b4cdbc06b4e2789d0246b121ce57935faea568c 100644
--- a/wetb/hawc2/hawc2_input_writer.py
+++ b/wetb/hawc2/hawc2_input_writer.py
@@ -1,61 +1,115 @@
+import itertools
+import jinja2
 import pandas as pd
-from wetb.hawc2.htc_file import HTCFile
-from wetb.hawc2.tests import test_files
+import click
 import os
-import itertools
-import importlib
+from pathlib import Path
 from pandas.core.base import PandasObject
-import numpy as np
-
+from wetb.hawc2.htc_file import HTCFile
 
-class HAWC2InputWriter():
 
+class HAWC2InputWriter(object):
+    """
+    Base HAWC2InputWriter object class, using the tagless system. Subclasses are:
+     - JinjaWriter
+    """
+    
     def __init__(self, base_htc_file, **kwargs):
         self.base_htc_file = base_htc_file
+        self.content = None
         for k, v in kwargs.items():
             setattr(self, k, v)
 
-    def from_pandas(self, dataFrame, write_input_files=True):
+
+    def __call__(self, out_fn, **kwargs):
+        return self.write(out_fn, **kwargs)
+
+
+    def from_pandas(self, dataFrame):
+        """
+        Loads a DataFrame of contents from a PandasObject
+        
+        Args:
+            dataFrame (Pandas.DataFrame of PandasObject): Dataframe of contents
+        
+        Returns:
+            self (HAWC2InputWriter): Instance of self to enable chaining.
+        """
+        
         if not isinstance(dataFrame, PandasObject) and hasattr(dataFrame, 'to_pandas'):
             dataFrame = dataFrame.to_pandas()
+            
         if isinstance(dataFrame, pd.Panel):
-            return {n: self.from_pandas(df.dropna(how='all'), write_input_files)for n, df in dataFrame.iteritems()}
+            return {n: self.from_pandas(df.dropna(how='all'))for n, df in dataFrame.iteritems()}
 
         else:
             assert 'Name' in dataFrame
-            if write_input_files:
-                for _, row in dataFrame.iterrows():
-                    self.write_input_files(**row.to_dict())
-        return dataFrame
-
-    def from_excel(self, excel_file, write_input_files=True):
-        return self.from_pandas(pd.read_excel(excel_file), write_input_files)
-
-    def from_CVF(self, constants, variables, functions, write_input_files=True):
-        attributes = list(constants.keys()) + list(variables.keys()) + list(functions.keys())
-        tags = []
+            
+        self.contents = dataFrame
+        
+        return self
+
+
+    def from_excel(self, excel_file):
+        """
+        Loads a DataFrame of contents from an excel file
+        
+        Args:
+            excel_file (str, Pathlib.Path): path to excel file.
+        
+        Returns:
+            self (HAWC2InputWriter): Instance of self to enable chaining.
+        """
+        
+        self.contents = pd.read_excel(excel_file)
+        return self
+
+
+    def from_CVF(self, constants, variables, functionals):
+        """
+        Produces a DataFrame of contents given dictionaries of constants, variables
+        and functionals. 
+
+        Args: 
+            constants (dict): content, value pairs for which the value remains
+                constant for all generated htc files. 
+            variables (dict of list): content, list pairs for which all combinations
+                of the listed values are generated in the htc files. 
+            functionals (dict of functions): content, function pairs for which the
+                content is dependent on the constants and variables contents.
+
+        Returns: self (HAWC2InputWriter): Instance of self to enable chaining.
+        """
+        
+        attributes = list(constants.keys()) + list(variables.keys()) + list(functionals.keys())
+        contents = []
+        
         for var in itertools.product(*list(variables.values())):
             this_dict = dict(constants)
             var_dict = dict(zip(variables.keys(), var))
             this_dict.update(var_dict)
-            for key, func in functions.items():
+            
+            for key, func in functionals.items():
                 this_dict[key] = func(this_dict)
 
-            tags.append(list(this_dict.values()))
+            contents.append(list(this_dict.values()))
 
-        df = pd.DataFrame(tags, columns=attributes)
-        return self.from_pandas(df, write_input_files)
+        self.contents = pd.DataFrame(contents, columns=attributes)
+        
+        return self
 
-    def from_definition(self, definition_file, write_input_files=True):
-        file = os.path.splitext(definition_file)[0]
-        module = os.path.relpath(file, os.getcwd()).replace(os.path.sep, ".")
-        def_mod = importlib.import_module(module)
-        return self.from_CVF(def_mod.constants, def_mod.variables, def_mod.functions, write_input_files)
 
-    def __call__(self, name, folder='', **kwargs):
-        return self.write_input_files(name, folder, **kwargs)
-
-    def write_input_files(self, Name, Folder='', DLC=None, **kwargs):
+    def write(self, out_fn, **kwargs):
+        ''' 
+        Renders a single htc file for a given set of 'tagless' contents. 
+        args: 
+        out_fn (str or pathlib.Path): The output filename where to render
+        the jinja template. 
+        params (pd.DataFrame or pd.Series) : The input contents.
+        '''
+        # if isinstance(params, PandasObject):
+        #     params = params.to_dict()
+            
         htc = HTCFile(self.base_htc_file)
         for k, v in kwargs.items():
             k = k.replace('/', '.')
@@ -63,42 +117,87 @@ class HAWC2InputWriter():
                 line = htc[k]
                 v = str(v).strip().replace(",", " ")
                 line.values = v.split()
+            
+            elif k in ['Name', 'Folder', 'DLC']:
+                continue
+            
             else:
                 getattr(self, 'set_%s' % k)(htc, **kwargs)
-        htc.set_name(Name, Folder)
-        htc.save()
-        args = {'Name': Name, 'Folder': Folder}
-        args.update(kwargs)
-
-        return pd.Series(args)
-
-
-def main():
-    if __name__ == '__main__':
-
-        path = os.path.dirname(test_files.__file__) + '/simulation_setup/DTU10MWRef6.0/'
-        base_htc_file = path + 'htc/DTU_10MW_RWT.htc'
-        h2writer = HAWC2InputWriter(base_htc_file)
-
-        # pandas2htc
-        df = pd.DataFrame({'wind.wsp': [4, 6, 8], 'Name': ['c1', 'c2', 'c3']})
-        h2writer.from_pandas(df)
-
-        # HAWC2InputWriter
-        class MyWriter(HAWC2InputWriter):
-            def set_time(self, htc, time, **_):
-                htc.set_time(self.time_start, self.time_start + time)
-
-        myWriter = MyWriter(base_htc_file, time_start=100)
-        myWriter("t1", 'tmp', time=600, **{"wind.wsp": 10})
-
-        # from_CVF (constants, variables and functions)
-        constants = {'simulation.time_stop': 100}
-        variables = {'wind.wsp': [4, 6, 8],
-                     'wind.tint': [0.1, 0.15, 0.2]}
-        functions = {'Name': lambda x: 'sim_wsp' + str(x['wind.wsp']) + '_ti' + str(x['wind.tint'])}
-        df = h2writer.from_CVF(constants, variables, functions, write_input_files=False)
-        print(df)
-
 
-main()
+            htc.save(out_fn)
+    
+        
+    def write_all(self, out_dir):
+        ''' 
+        Renders all htc files for the set of contents.
+        args:
+            out_dir (str or pathlib.Path): The directory where the htc files are generated.
+        '''
+        out_dir = Path(out_dir)
+        out_dir.mkdir(parents=True, exist_ok=True)
+        
+        N = len(self.contents)
+        
+        print(f'Generating {N} htc files in directory: {out_dir}')
+        
+        with click.progressbar(self.contents.iterrows(), length=N) as bar:
+            for _, row in bar:
+                self.write(out_dir/'{}.htc'.format(row.Name), **row.to_dict())
+
+                
+                
+class JinjaWriter(HAWC2InputWriter):
+    """
+    Subclass of the HAWC2InputWriter object. Generates htc files using contents compatible with jinja2.
+    """
+    def write(self, out_fn, **kwargs):
+        params = pd.Series(kwargs)
+        with open(self.base_htc_file) as f:
+            template = jinja2.Template(f.read())
+            
+            out = template.render(params)
+            
+            with open(out_fn, 'w') as f:
+                f.write(out)
+
+
+
+
+if __name__ == '__main__':
+    from wetb.hawc2.tests import test_files
+    
+    ### HAWC2InputWriter example
+    path = os.path.dirname(test_files.__file__) + '/simulation_setup/DTU10MWRef6.0/'
+    base_htc_file = path + 'htc/DTU_10MW_RWT.htc'
+
+    class MyWriter(HAWC2InputWriter):
+        def set_time(self, htc, time, **_):
+            htc.set_time(self.time_start, self.time_start + time)
+
+    myWriter = MyWriter(base_htc_file, time_start=100)
+    myWriter('tmp/t1.htc', Name="t1", time=600, **{"wind.wsp": 10})
+
+    
+    
+    ### JinjaWriter example with constants, variables, and functionals
+    Constants = {
+        'time_stop': 100,
+    }
+
+    Variables = {
+        'wsp'  : [4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24],
+        'tint' : [0.05, 0.10, 0.15],
+
+    }
+
+    Functionals = {
+        'Name': lambda x: 'dtu10mw_wsp' + str(x['wsp']) + '_ti' + str(x['tint']),
+    }
+    
+    htc_template_fn = os.path.dirname(test_files.__file__) + '/simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc.j2'
+    writer = JinjaWriter(htc_template_fn)
+    writer.from_CVF(Constants, Variables, Functionals)
+    print(writer.contents)
+    
+    # write all htc files to folder 'htc'
+    writer.write_all('htc')
diff --git a/wetb/hawc2/tests/test_files/simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc.j2 b/wetb/hawc2/tests/test_files/simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc.j2
new file mode 100644
index 0000000000000000000000000000000000000000..5f08607702198ecfed5723b281fd82db4278e730
--- /dev/null
+++ b/wetb/hawc2/tests/test_files/simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc.j2
@@ -0,0 +1,652 @@
+;DTU_10MW_RWT, version 5, 04-21-15, anyd
+;
+  begin simulation;
+    time_stop	{{time_stop}};
+    solvertype	1;		(newmark)
+    on_no_convergence	continue;
+    convergence_limits	1000 1 1e-07;
+    logfile	./log/{{Name}}.log;
+    begin newmark;
+      deltat	0.02;
+    end newmark;
+  end simulation;		;
+;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+  begin new_htc_structure;		;   beam_output_file_name  ./log/DTU_10MW_RWT_beam.dat;                    Optional - Calculated beam properties of the bodies are written to file
+;   body_output_file_name  ./log/DTU_10MW_RWT_body.dat;                    Optional - Body initial position and orientation are written to file
+;   body_eigenanalysis_file_name ./eig/DTU_10MW_RWT_body_eigen.dat;
+;   structure_eigenanalysis_file_name ./eig/DTU_10MW_RWT_strc_eigen.dat ;
+;-------------------------------------------------------------------------------------------------------------------------------
+;-------------------------------------------------------------------------------------------------------------------------------
+    begin main_body;		tower 115m
+      name	tower;
+      type	timoschenko;
+      nbodies	1;
+      node_distribution	c2_def;
+      damping_posdef	0 0 0 0.00412 0.00412 0.00045;		Mx My Mz Kx Ky Kz , M´s raises overall level, K´s raises high freguency level "tuned by Larh"
+      begin timoschenko_input;
+        filename	./data/dtu_10mw_rwt_tower_st.dat;
+        set	1 1;
+      end timoschenko_input;
+      begin c2_def;		Definition of centerline (main_body coordinates)
+        nsec	11;
+        sec	1 0 0 0 0;		x,y,z,twist
+        sec	2 0 0 -11.5 0;
+        sec	3 0 0 -23 0;
+        sec	4 0 0 -34.5 0;
+        sec	5 0 0 -46 0;
+        sec	6 0 0 -57.5 0;
+        sec	7 0 0 -69 0;
+        sec	8 0 0 -80.5 0;
+        sec	9 0 0 -92 0;
+        sec	10 0 0 -103.5 0;
+        sec	11 0 0 -115.63 0;
+      end c2_def;
+    end main_body;		;
+    begin main_body;
+      name	towertop;
+      type	timoschenko;
+      nbodies	1;
+      node_distribution	c2_def;
+      damping_posdef	0 0 0 0.007 0.007 0.007;		"changed by Larh"
+      concentrated_mass	2 0 2.687 0.30061 446040 4106000 410600 4106000;		Nacelle mass and inertia "corrected by Anyd 25/4/13"
+      begin timoschenko_input;
+        filename	./data/dtu_10mw_rwt_towertop_st.dat;
+        set	1 2;
+      end timoschenko_input;
+      begin c2_def;		Definition of centerline (main_body coordinates)
+        nsec	2;
+        sec	1 0 0 0 0;		x,y,z,twist
+        sec	2 0 0 -2.75 0;
+      end c2_def;
+    end main_body;		;
+    begin main_body;
+      name	shaft;
+      type	timoschenko;
+      nbodies	1;
+      node_distribution	c2_def;
+      damping_posdef	0 0 0 0.000465 0.000465 0.003983;		"tuned by Anyd 23/5/13 to 31.45 log decr. damping for free free with stiff rotor and tower"
+      concentrated_mass	1 0 0 0 0 0 0 3751000;		generator equivalent slow shaft "re_tuned by Anyd 20/2/13"
+      concentrated_mass	5 0 0 0 105520 0 0 325700;		hub mass and inertia;	"re_tuned by Anyd 20/2/13"
+      begin timoschenko_input;
+        filename	./data/dtu_10mw_rwt_shaft_st.dat;
+        set	1 1;
+      end timoschenko_input;
+      begin c2_def;		Definition of centerline (main_body coordinates)
+        nsec	5;
+        sec	1 0 0 0 0;		Tower top x,y,z,twist
+        sec	2 0 0 1.5 0;
+        sec	3 0 0 3 0;
+        sec	4 0 0 4.4 0;		Main bearing
+        sec	5 0 0 7.1 0;		Rotor centre
+      end c2_def;
+    end main_body;		;
+    begin main_body;
+      name	hub1;
+      type	timoschenko;
+      nbodies	1;
+      node_distribution	c2_def;
+      damping_posdef	0 0 0 3e-06 3e-06 2e-05;		"changed by Larh"
+      begin timoschenko_input;
+        filename	./data/dtu_10mw_rwt_hub_st.dat;
+        set	1 2;
+      end timoschenko_input;
+      begin c2_def;		Definition of centerline (main_body coordinates)
+        nsec	2;
+        sec	1 0 0 0 0;		x,y,z,twist
+        sec	2 0 0 2.8 0;
+      end c2_def;
+    end main_body;		;
+    begin main_body;
+      name	hub2;
+      copy_main_body	hub1;
+    end main_body;		;
+    begin main_body;
+      name	hub3;
+      copy_main_body	hub1;
+    end main_body;		;
+    begin main_body;
+      name	blade1;
+      type	timoschenko;
+      nbodies	10;
+      node_distribution	c2_def;
+      damping_posdef	0 0 0 0.00153 0.00255 0.00033;		" 3% damping tuned by tkim 23/03/13 unable to fit 3rd and higher mode"
+      begin timoschenko_input;
+        filename	./data/dtu_10mw_rwt_blade_st.dat;
+        set	1 1;		set subset
+      end timoschenko_input;
+      begin c2_def;		Definition of centerline (main_body coordinates)
+        nsec	27;
+        sec	1 0 7.006e-05 4.44089e-16 -14.5;
+        sec	2 -2.06477e-05 -0.0122119 3 -14.5;
+        sec	3 -0.0072881 -0.0249251 6 -14.4851;
+        sec	4 -0.0189235 -0.0273351 7.00004 -14.461;
+        sec	5 -0.0541282 -0.0282163 8.70051 -14.3388;
+        sec	6 -0.126633 -0.021321 10.402 -14.0201;
+        sec	7 -0.225666 -0.0128378 12.2046 -13.3904;
+        sec	8 -0.288563 -0.00770659 13.2065 -12.9371;
+        sec	9 -0.399194 -0.00488317 15.01 -11.9445;
+        sec	10 -0.576634 -0.0180296 18.2151 -9.98243;
+        sec	11 -0.707136 -0.0501772 21.4178 -8.45147;
+        sec	12 -0.791081 -0.0941228 24.6189 -7.46417;
+        sec	13 -0.837195 -0.14888 27.8193 -6.72916;
+        sec	14 -0.853948 -0.214514 31.0194 -6.08842;
+        sec	15 -0.849367 -0.290618 34.2197 -5.49322;
+        sec	16 -0.79392 -0.462574 40.2204 -4.39222;
+        sec	17 -0.716284 -0.688437 46.6217 -3.09315;
+        sec	18 -0.634358 -0.960017 53.0232 -1.75629;
+        sec	19 -0.553179 -1.28424 59.4245 -0.50065;
+        sec	20 -0.475422 -1.66402 65.8255 0.601964;
+        sec	21 -0.40318 -2.10743 72.2261 1.5556;
+        sec	22 -0.330085 -2.6563 79.0266 2.51935;
+        sec	23 -0.31014 -2.78882 80.5267 2.7295;
+        sec	24 -0.286719 -2.92517 82.0271 2.93201;
+        sec	25 -0.255823 -3.06577 83.5274 3.11874;
+        sec	26 -0.207891 -3.20952 85.0277 3.28847;
+        sec	27 -0.089894 -3.33685 86.3655 3.42796;
+      end c2_def;
+    end main_body;		;
+    begin main_body;
+      name	blade2;
+      copy_main_body	blade1;
+    end main_body;		;
+    begin main_body;
+      name	blade3;
+      copy_main_body	blade1;
+    end main_body;		;-------------------------------------------------------------------------------------------------------------------------------
+;
+    begin orientation;
+      begin base;
+        body	tower;
+        inipos	0 0 0;		initial position of node 1
+        body_eulerang	0 0 0;
+      end base;		;
+      begin relative;
+        body1	tower last;
+        body2	towertop 1;
+        body2_eulerang	0 0 0;
+      end relative;		;
+      begin relative;
+        body1	towertop last;
+        body2	shaft 1;
+        body2_eulerang	90 0 0;
+        body2_eulerang	5 0 0;		5 deg tilt angle
+        body2_eulerang	0 0 0;
+        mbdy2_ini_rotvec_d1	0 0 -1 0.7;		mbdy2_ini_rotvec_d1 0.0 0.0 -1.0 [init_wr];
+      end relative;		;
+      begin relative;
+        body1	shaft last;
+        body2	hub1 1;
+        body2_eulerang	-90 0 0;
+        body2_eulerang	0 180 0;
+        body2_eulerang	2.5 0 0;		2.5deg cone angle
+      end relative;		;
+      begin relative;
+        body1	shaft last;
+        body2	hub2 1;
+        body2_eulerang	-90 0 0;
+        body2_eulerang	0 60 0;
+        body2_eulerang	2.5 0 0;		2.5deg cone angle
+      end relative;		;
+      begin relative;
+        body1	shaft last;
+        body2	hub3 1;
+        body2_eulerang	-90 0 0;
+        body2_eulerang	0 -60 0;
+        body2_eulerang	2.5 0 0;		2.5deg cone angle
+      end relative;		;
+      begin relative;
+        body1	hub1 last;
+        body2	blade1 1;
+        body2_eulerang	0 0 0;
+      end relative;		;
+      begin relative;
+        body1	hub2 last;
+        body2	blade2 1;
+        body2_eulerang	0 0 0;
+      end relative;		;
+      begin relative;
+        body1	hub3 last;
+        body2	blade3 1;
+        body2_eulerang	0 0 0;
+      end relative;		;
+    end orientation;		;-------------------------------------------------------------------------------------------------------------------------------
+    begin constraint;		;
+      begin fix0;		fixed to ground in translation and rotation of node 1
+        body	tower;
+      end fix0;		;
+      begin fix1;
+        body1	tower last;
+        body2	towertop 1;
+      end fix1;		;
+      begin bearing1;		free bearing
+        name	shaft_rot;
+        body1	towertop last;
+        body2	shaft 1;
+        bearing_vector	2 0 0 -1;		x=coo (0=global.1=body1.2=body2) vector in body2 coordinates where the free rotation is present
+      end bearing1;		;
+      begin fix1;
+        body1	shaft last;
+        body2	hub1 1;
+      end fix1;		;
+      begin fix1;
+        body1	shaft last;
+        body2	hub2 1;
+      end fix1;		;
+      begin fix1;
+        body1	shaft last;
+        body2	hub3 1;
+      end fix1;		;
+      begin bearing2;
+        name	pitch1;
+        body1	hub1 last;
+        body2	blade1 1;
+        bearing_vector	2 0 0 -1;
+      end bearing2;		;
+      begin bearing2;
+        name	pitch2;
+        body1	hub2 last;
+        body2	blade2 1;
+        bearing_vector	2 0 0 -1;
+      end bearing2;		;
+      begin bearing2;
+        name	pitch3;
+        body1	hub3 last;
+        body2	blade3 1;
+        bearing_vector	2 0 0 -1;
+      end bearing2;
+    end constraint;		;
+  end new_htc_structure;		;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+  begin wind;
+    density	1.225;
+    wsp	{{wsp}};
+    tint	{{tint}};
+    horizontal_input	1;
+    windfield_rotations	0 0 0;		yaw, tilt, rotation
+    center_pos0	0 0 -119;		hub heigth
+    shear_format	1 0.2;
+    turb_format	1;		0=none, 1=mann,2=flex
+    tower_shadow_method	0;		0=none, 1=potential flow, 2=jet
+    scale_time_start	0;
+    wind_ramp_factor	0 40 0.6 1;		; Steps ;
+    wind_ramp_abs	140 141 0 1;		wsp. after the step:  5.0
+    wind_ramp_abs	181 182 0 1;		wsp. after the step:  6.0
+    wind_ramp_abs	222 223 0 1;		wsp. after the step:  7.0
+    wind_ramp_abs	263 264 0 1;		wsp. after the step:  8.0
+    wind_ramp_abs	304 305 0 1;		wsp. after the step:  9.0
+    wind_ramp_abs	345 346 0 1;		wsp. after the step: 10.0
+    wind_ramp_abs	386 387 0 1;		wsp. after the step: 11.0
+    wind_ramp_abs	427 428 0 1;		wsp. after the step: 12.0
+    wind_ramp_abs	468 469 0 1;		wsp. after the step: 13.0
+    wind_ramp_abs	509 510 0 1;		wsp. after the step: 14.0
+    wind_ramp_abs	550 551 0 1;		wsp. after the step: 15.0
+    wind_ramp_abs	591 592 0 1;		wsp. after the step: 16.0
+    wind_ramp_abs	632 633 0 1;		wsp. after the step: 17.0
+    wind_ramp_abs	673 674 0 1;		wsp. after the step: 18.0
+    wind_ramp_abs	714 715 0 1;		wsp. after the step: 19.0
+    wind_ramp_abs	755 756 0 1;		wsp. after the step: 20.0
+    wind_ramp_abs	796 797 0 1;		wsp. after the step: 21.0
+    wind_ramp_abs	837 838 0 1;		wsp. after the step: 22.0
+    wind_ramp_abs	878 879 0 1;		wsp. after the step: 23.0
+    wind_ramp_abs	919 920 0 1;		wsp. after the step: 24.0
+    wind_ramp_abs	960 961 0 1;		wsp. after the step: 25.0
+;
+    begin tower_shadow_potential_2;
+      tower_mbdy_link	tower;
+      nsec	2;
+      radius	0 4.15;
+      radius	115.63 2.75;
+    end tower_shadow_potential_2;
+  end wind;		;
+  begin aerodrag;
+    begin aerodrag_element;
+      mbdy_name	tower;
+      aerodrag_sections	uniform 10;
+      nsec	2;
+      sec	0 0.6 8.3;		tower bottom
+      sec	115.63 0.6 5.5;		tower top
+    end aerodrag_element;		;
+    begin aerodrag_element;		Nacelle drag side
+      mbdy_name	shaft;
+      aerodrag_sections	uniform 2;
+      nsec	2;
+      sec	0 0.8 10;
+      sec	7.01 0.8 10;
+    end aerodrag_element;
+  end aerodrag;		;
+  begin aero;
+    nblades	3;
+    hub_vec	shaft -3;		rotor rotation vector (normally shaft composant directed from pressure to sustion side)
+    link	1 mbdy_c2_def blade1;
+    link	2 mbdy_c2_def blade2;
+    link	3 mbdy_c2_def blade3;
+    ae_filename	./data/dtu_10mw_rwt_ae.dat;
+    pc_filename	./data/dtu_10mw_rwt_pc.dat;
+    induction_method	1;		0=none, 1=normal
+    aerocalc_method	1;		0=ingen aerodynamic, 1=med aerodynamic
+    aerosections	50;		def. 50
+    ae_sets	1 1 1;
+    tiploss_method	1;		0=none, 1=prandtl
+    dynstall_method	2;		0=none, 1=stig øye method,2=mhh method
+;
+  end aero;		;-------------------------------------------------------------------------------------------------
+  begin dll;		;
+    begin type2_dll;
+      name	risoe_controller;
+      filename	./control/dtu_we_controller.dll;
+      dll_subroutine_init	init_regulation;
+      dll_subroutine_update	update_regulation;
+      arraysizes_init	52 1;
+      arraysizes_update	12 100;
+      begin init;		; Overall parameters
+        constant	1 10000;		Rated power [kW]
+        constant	2 0.5236;		Minimum rotor speed [rad/s]
+        constant	3 1.005;		Rated rotor speed [rad/s]
+        constant	4 15600000;		Maximum allowable generator torque [Nm]
+        constant	5 100;		Minimum pitch angle, theta_min [deg],
+; if |theta_min|>90, then a table of <wsp,theta_min> is read ;
+; from a file named 'wptable.n', where n=int(theta_min)
+        constant	6 82;		Maximum pitch angle [deg]
+        constant	7 10;		Maximum pitch velocity operation [deg/s]
+        constant	8 0.4;		Frequency of generator speed filter [Hz]
+        constant	9 0.7;		Damping ratio of speed filter [-]
+        constant	10 1.92;		Frequency of free-free DT torsion mode [Hz], if zero no notch filter used
+; Partial load control parameters
+        constant	11 11750000;		Optimal Cp tracking K factor [Nm/(rad/s)^2], ;
+; Qg=K*Omega^2, K=eta*0.5*rho*A*Cp_opt*R^3/lambda_opt^3
+        constant	12 70840000;		Proportional gain of torque controller [Nm/(rad/s)]
+        constant	13 15900000;		Integral gain of torque controller [Nm/rad]
+        constant	14 0;		Differential gain of torque controller [Nm/(rad/s^2)]
+;     Full load control parameters
+        constant	15 2;		Generator control switch [1=constant power, 2=constant torque]
+        constant	16 1.304;		Proportional gain of pitch controller [rad/(rad/s)]
+        constant	17 0.3511;		Integral gain of pitch controller [rad/rad]
+        constant	18 0;		Differential gain of pitch controller [rad/(rad/s^2)]
+        constant	19 4e-09;		Proportional power error gain [rad/W]
+        constant	20 4e-09;		Integral power error gain [rad/(Ws)]
+        constant	21 11.35;		Coefficient of linear term in aerodynamic gain scheduling, KK1 [deg]
+        constant	22 400.7;		Coefficient of quadratic term in aerodynamic gain scheduling, KK2 [deg^2] &
+; (if zero, KK1 = pitch angle at double gain)
+        constant	23 1.3;		Relative speed for double nonlinear gain [-]
+;     Cut-in simulation parameters
+        constant	24 -1;		Cut-in time [s]
+        constant	25 1;		Time delay for soft start of torque [1/1P]
+;     Cut-out simulation parameters
+        constant	26 1500;		Cut-out time [s]
+        constant	27 5;		Time constant for linear torque cut-out [s]
+        constant	28 1;		Stop type [1=normal, 2=emergency]
+        constant	29 1;		Time delay for pitch stop after shut-down signal [s]
+        constant	30 3;		Maximum pitch velocity during initial period of stop [deg/s]
+        constant	31 3;		Time period of initial pitch stop phase [s] (maintains pitch speed specified in constant 30)
+        constant	32 4;		Maximum pitch velocity during final phase of stop [deg/s]
+;     Expert parameters (keep default values unless otherwise given)
+        constant	33 2;		Lower angle above lowest minimum pitch angle for switch [deg]
+        constant	34 2;		Upper angle above lowest minimum pitch angle for switch [deg], if equal then hard switch
+        constant	35 95;		Ratio between filtered speed and reference speed for fully open torque limits [%]
+        constant	36 2;		Time constant of 1st order filter on wind speed used for minimum pitch [1/1P]
+        constant	37 1;		Time constant of 1st order filter on pitch angle used for gain scheduling [1/1P]
+;     Drivetrain damper
+        constant	38 0;		Proportional gain of active DT damper [Nm/(rad/s)], requires frequency in input 10
+;	  Over speed
+        constant	39 25;		Overspeed percentage before initiating turbine controller alarm (shut-down) [%]
+;     Additional non-linear pitch control term (not used when all zero)
+        constant	40 0;		Err0 [rad/s]
+        constant	41 0;		ErrDot0 [rad/s^2]
+        constant	42 0;		PitNonLin1 [rad/s]
+;     Storm control command
+        constant	43 28;		Wind speed 'Vstorm' above which derating of rotor speed is used [m/s]
+        constant	44 28;		Cut-out wind speed (only used for derating of rotor speed in storm) [m/s]
+;     Safety system parameters
+        constant	45 300;		Overspeed percentage before initiating safety system alarm (shut-down) [%]
+        constant	46 1.5;		Max low-pass filtered tower top acceleration level [m/s^2]
+;     Turbine parameter
+        constant	47 178;		Nominal rotor diameter [m]
+;     Parameters for rotor inertia reduction in variable speed region
+        constant	48 0;		Proportional gain on rotor acceleration in variable speed region [Nm/(rad/s^2)] (not used when zero)
+;     Parameters for alternative partial load controller with PI regulated TSR tracking
+        constant	49 0;		Optimal tip speed ratio [-] (only used when K=constant 11 = 0 otherwise  Qg=K*Omega^2 is used)
+;     Parameters for adding aerodynamic drivetrain damping on gain scheduling
+        constant	50 0;		Proportional gain of aerodynamic DT damping [Nm/(rad/s)]
+        constant	51 0;		Coefficient of linear term in aerodynamic DT damping scheduling, KK1 [deg]
+        constant	52 0;		Coefficient of quadratic term in aerodynamic DT damping scheduling, KK2 [deg^2]
+      end init;		;
+      begin output;
+        general time;						 [s]
+        constraint bearing1	shaft_rot 1 only 2;						 Drivetrain speed [rad/s]
+        constraint bearing2	pitch1 1 only 1;						 [rad]
+        constraint bearing2	pitch2 1 only 1;						 [rad]
+        constraint bearing2	pitch3 1 only 1;						 [rad]
+        wind free_wind	1 0 0 -119;						 Global coordinates at hub height
+        dll inpvec	2 2;						 Elec. power from generator servo .dll
+        dll inpvec	2 8;						 Grid state flag from generator servo .dll
+        mbdy state	acc tower 10 1 global only 1;						 Tower top x-acceleration [m/s^2]
+        mbdy state	acc tower 10 1 global only 2;						 Tower top y-acceleration [m/s^2]
+      end output;
+    end type2_dll;		;
+    begin type2_dll;
+      name	generator_servo;
+      filename	./control/generator_servo.dll;
+      dll_subroutine_init	init_generator_servo;
+      dll_subroutine_update	update_generator_servo;
+      arraysizes_init	7 1;
+      arraysizes_update	4 8;
+      begin init;
+        constant	1 20;		Frequency of 2nd order servo model of generator-converter system [Hz]
+        constant	2 0.9;		Damping ratio 2nd order servo model of generator-converter system [-]
+        constant	3 15600000;		Maximum allowable LSS torque (pull-out torque) [Nm]
+        constant	4 0.94;		Generator efficiency [-]
+        constant	5 1;		Gearratio [-]
+        constant	6 0;		Time for half value in softstart of torque [s]
+        constant	7 1500;		Time for grid loss [s]
+      end init;		;
+      begin output;
+        general time;						 Time [s]
+        dll inpvec	1 1;						 Electrical torque reference [Nm]
+        constraint bearing1	shaft_rot 1 only 2;						 Generator LSS speed [rad/s]
+        mbdy momentvec	shaft 1 1 shaft only 3;						 Shaft moment [kNm] (Qshaft)
+      end output;	
+;
+      begin actions;
+        mbdy	moment_int shaft 1 -3 shaft towertop 2;		Generator LSS torque [Nm]
+      end actions;
+    end type2_dll;		;
+    begin type2_dll;
+      name	mech_brake;
+      filename	./control/mech_brake.dll;
+      dll_subroutine_init	init_mech_brake;
+      dll_subroutine_update	update_mech_brake;
+      arraysizes_init	7 1;
+      arraysizes_update	4 6;
+      begin init;
+        constant	1 2727252;		Fully deployed maximum brake torque [Nm]
+        constant	2 100;		Parameter alpha used in Q = tanh(omega*alpha), typically 1e2/Omega_nom
+        constant	3 0.625;		Delay time for before brake starts to deploy [s] - from 5MW*1P_5/1P_10
+        constant	4 0.75;		Time for brake to become fully deployed [s]
+      end init;		;
+      begin output;
+        general time;						 Time [s]
+        constraint bearing1	shaft_rot 1 only 2;						 Generator LSS speed [rad/s]
+        dll inpvec	1 25;						 Command to deploy mechanical disc brake [0,1]
+      end output;	
+;
+      begin actions;
+        mbdy	moment_int shaft 1 3 shaft towertop 2;		Brake LSS torque [Nm]
+      end actions;
+    end type2_dll;		;
+    begin type2_dll;
+      name	servo_with_limits;
+      filename	./control/servo_with_limits.dll;
+      dll_subroutine_init	init_servo_with_limits;
+      dll_subroutine_update	update_servo_with_limits;
+      arraysizes_init	10 1;
+      arraysizes_update	5 9;
+      begin init;
+        constant	1 3;		Number of blades [-]
+        constant	2 1;		Frequency of 2nd order servo model of pitch system [Hz]
+        constant	3 0.7;		Damping ratio 2nd order servo model of pitch system [-]
+        constant	4 10;		Max. pitch speed [deg/s]
+        constant	5 15;		Max. pitch acceleration [deg/s^2]
+        constant	6 -5;		Min. pitch angle [deg]
+        constant	7 90;		Max. pitch angle [deg]
+        constant	8 1500;		Time for pitch runaway [s]
+        constant	9 1500;		Time for stuck blade 1 [s]
+        constant	10 0;		Angle of stuck blade 1 [deg]
+      end init;
+      begin output;
+        general time;						  Time                         [s]
+        dll inpvec	1 2;						  Pitch1 demand angle          [rad]
+        dll inpvec	1 3;						  Pitch2 demand angle          [rad]
+        dll inpvec	1 4;						  Pitch3 demand angle          [rad]
+        dll inpvec	1 26;						  Flag for emergency pitch stop         [0=off/1=on]
+      end output;	
+;
+      begin actions;
+        constraint	bearing2 angle pitch1;		Angle pitch1 bearing    [rad]
+        constraint	bearing2 angle pitch2;		Angle pitch2 bearing    [rad]
+        constraint	bearing2 angle pitch3;		Angle pitch3 bearing    [rad]
+      end actions;
+    end type2_dll;		;
+;	--- DLL for tower-blade tip distance -- ;
+    begin type2_dll;
+      name	disttowtip;
+      filename	./control/towclearsens.dll;
+      dll_subroutine_init	initialize;
+      dll_subroutine_update	update;
+      arraysizes_init	1 1;
+      arraysizes_update	12 4;
+      begin init;
+        constant	1 2.66;		Tower radius close to downward blade tip [m]
+      end init;
+      begin output;
+        mbdy state	pos tower 3 0.62 global;						 [1,2,3]. Tower position: 30.18 m
+        mbdy state	pos blade1 26 1 global;						 [4,5,6]
+        mbdy state	pos blade2 26 1 global;						 [7,8,9]
+        mbdy state	pos blade3 26 1 global;						 [10,11,12]
+      end output;
+    end type2_dll;
+  end dll;		;----------------------------------------------------------------------------------------------------------------------------------------------------------------
+;
+  begin output;
+    filename	./res/{{Name}};		; time 99.0 1000.0 ;
+; data_format  hawc_ascii;
+    data_format	hawc_binary;
+    buffer	1;		;
+    general time;
+    constraint bearing1	shaft_rot 2;			 angle and angle velocity
+    constraint bearing2	pitch1 5;			    angle and angle velocity
+    constraint bearing2	pitch2 5;			    angle and angle velocity
+    constraint bearing2	pitch3 5;			    angle and angle velocity
+    aero omega;
+    aero torque;
+    aero power;
+    aero thrust;
+    wind free_wind	1 0 0 -119;			 local wind at fixed position: coo (1=global,2=non-rotation rotor coo.), pos x, pos y, pos z
+; Moments:
+    mbdy momentvec	tower 1 1 tower # tower base;
+    mbdy momentvec	tower 10 2 tower # tower yaw bearing;
+    mbdy momentvec	shaft 4 1 shaft # main bearing;
+    mbdy momentvec	blade1 2 2 blade1 # blade 1 root;
+    mbdy momentvec	blade2 2 2 blade2 # blade 2 root;
+    mbdy momentvec	blade3 2 2 blade3 # blade 3 root;
+    mbdy momentvec	blade1 13 1 local # blade 1 50% local e coo;
+    mbdy momentvec	blade2 13 1 local # blade 2 50% local e coo;
+    mbdy momentvec	blade3 13 1 local # blade 3 50% local e coo;	
+; Displacements and accellerations
+    mbdy state	pos tower 10 1 global only 1 # tower top fa displ;
+    mbdy state	pos tower 10 1 global only 2 # tower top ss displ;
+    mbdy state	acc tower 10 1 global only 1 # tower top fa acc;
+    mbdy state	acc tower 10 1 global only 2 # tower top ss acc;	
+;
+    mbdy state	pos blade1 26 1 blade1 # blade 1 tip pos;
+    mbdy state	pos blade2 26 1 blade2 # blade 2 tip pos;
+    mbdy state	pos blade3 26 1 blade3 # blade 3 tip pos;
+    mbdy state	pos blade1 26 1 global # gl blade 1 tip pos;	
+; - Monitor Aerodynamics - ;
+    aero windspeed	3 1 1 72.5;
+    aero alfa	1 72.5;
+    aero alfa	2 72.5;
+    aero alfa	3 72.5;
+    aero cl	1 72.5;
+    aero cl	2 72.5;
+    aero cl	3 72.5;
+    aero cd	1 72.5;
+    aero cd	2 72.5;
+    aero cd	3 72.5;	
+; - Main Controller -
+; Output to controller
+; dll outvec 1 1 # time;
+; dll outvec 1 2 # slow speed shaft rad/s;
+; dll outvec 1 3 # pitch angle 1;
+; dll outvec 1 4 # pitch angle 2;
+; dll outvec 1 5 # pitch angle 3;
+; dll outvec 1 6 # WSP_x_global;
+; dll outvec 1 7 # WSP_y_global;
+; dll outvec 1 8 # WSP_z_global;
+; dll outvec 1 9 # Elec. pwr ;
+; dll outvec 1 10 # Grid flag ;
+; Input from controller
+    dll inpvec	1 1 # generator torque reference [nm];
+    dll inpvec	1 2 # pitch angle reference of blade 1 [rad];
+    dll inpvec	1 3 # pitch angle reference of blade 2 [rad];
+    dll inpvec	1 4 # pitch angle reference of blade 3 [rad];
+    dll inpvec	1 5 # power reference [w];
+    dll inpvec	1 6 # filtered wind speed [m/s];
+    dll inpvec	1 7 # filtered rotor speed [rad/s];
+    dll inpvec	1 8 # filtered rotor speed error for torque [rad/s];
+    dll inpvec	1 9 # bandpass filtered rotor speed [rad/s];
+    dll inpvec	1 10 # proportional term of torque contr. [nm];
+    dll inpvec	1 11 # integral term of torque controller [nm];
+    dll inpvec	1 12 # minimum limit of torque [nm];
+    dll inpvec	1 13 # maximum limit of torque [nm];
+    dll inpvec	1 14 # torque limit switch based on pitch [-];
+    dll inpvec	1 15 # filtered rotor speed error for pitch [rad/s];
+    dll inpvec	1 16 # power error for pitch [w];
+    dll inpvec	1 17 # proportional term of pitch controller [rad];
+    dll inpvec	1 18 # integral term of pitch controller [rad];
+    dll inpvec	1 19 # minimum limit of pitch [rad];
+    dll inpvec	1 20 # maximum limit of pitch [rad];
+    dll inpvec	1 21 # torque reference from dt dammper [nm];
+    dll inpvec	1 22 # status signal [-];
+    dll inpvec	1 23 # total added pitch rate [rad/s];
+    dll inpvec	1 24 # filtered mean pitch for gain sch [rad];
+    dll inpvec	1 25 # flag for mechnical brake [0=off/1=on];
+    dll inpvec	1 26 # flag for emergency pitch stop [0=off/1=on];
+    dll inpvec	1 27 # lp filtered acceleration level [m/s^2];	
+; ; Output to generator model
+; dll outvec 2 1  # time ;
+; dll outvec 2 2  # Electrical torque reference [Nm] ;
+; dll outvec 2 3  # omega LSS ;
+; Input from generator model
+    dll inpvec	2 1 # mgen lss [nm];
+    dll inpvec	2 2 # pelec w;
+    dll inpvec	2 3 # mframe;
+    dll inpvec	2 4 # mgen hss;
+    dll inpvec	2 5 # generator pmech kw;
+    dll inpvec	2 6 # filtered gen speed;
+    dll inpvec	2 7 # elec. pwr;
+    dll inpvec	2 8 # grid flag;	
+; Output to mechanical brake
+    dll inpvec	3 1 # brake torque [nm];	
+; ; Input from mechanical brake
+; dll outvec 3 1 # Time [s] ;
+; dll outvec 3 2 # Generator LSS speed [rad/s] ;
+; dll outvec 3 3 # Deploy brake ;
+; ; Output to pitch servo
+; dll outvec 4 1 # time;
+; dll outvec 4 2 # pitchref 1;
+; dll outvec 4 3 # pitchref 2;
+; dll outvec 4 4 # pitchref 3;
+; dll outvec 4 5 # Emerg. stop;
+; Input from pitch servo
+    dll inpvec	4 1 # pitch 1;
+    dll inpvec	4 2 # pitch 2;
+    dll inpvec	4 3 # pitch 3;	
+; Check tower clearence
+    dll inpvec	5 1 # bltip tow min d [m];
+  end output;	
+;
+  begin output_at_time aero 0.1;
+    filename	res/at;
+    twist 1;
+    chord 1;
+  end output_at_time;
+  exit;
diff --git a/wetb/hawc2/tests/test_hawc2_input_writer.py b/wetb/hawc2/tests/test_hawc2_input_writer.py
index da9fbceef6b346c0eb841c45d9e0cd024c1b89bb..b81657d2fcc0c900d42555b4496edbc8e54be8f8 100644
--- a/wetb/hawc2/tests/test_hawc2_input_writer.py
+++ b/wetb/hawc2/tests/test_hawc2_input_writer.py
@@ -34,6 +34,7 @@ def test_pandas2htc(h2writer):
     wsp_lst = [4, 6, 8]
     df = pd.DataFrame({'wind.wsp': wsp_lst, 'Name': ['c1', 'c2', 'c3'], 'Folder': ['tmp', 'tmp', 'tmp']})
     h2writer.from_pandas(df)
+    h2writer.write_all(path + 'htc/tmp/')
     for i, wsp in enumerate(wsp_lst, 1):
         htc = HTCFile(path + "htc/tmp/c%d.htc" % i)
         assert htc.wind.wsp[0] == wsp
@@ -41,6 +42,7 @@ def test_pandas2htc(h2writer):
 
 def test_excel2htc(h2writer):
     h2writer.from_excel(os.path.dirname(test_files.__file__) + "/htc_input_table.xlsx")
+    h2writer.write_all(path + 'htc/tmp/')
     for i, wsp in enumerate([4, 6], 1):
         htc = HTCFile(path + "htc/tmp/a%d.htc" % i)
         assert htc.wind.wsp[0] == wsp
@@ -60,11 +62,11 @@ def test_hawc2_writer():
             htc.set_time(self.time_start, self.time_start + time)
 
     myWriter = MyWriter(htc_base_file, time_start=100)
-    ps = myWriter("w1", 'tmp', time=600, **{"wind.wsp": 10})
+    ps = myWriter.write(path + "htc/tmp/w1.htc", **{"Name":"w1", "wind.wsp": 10, "time":600})
     htc = HTCFile(path + "htc/tmp/w1.htc")
     assert htc.simulation.time_stop[0] == 700
     assert htc.wind.wsp[0] == 10
-    assert (ps.Name == 'w1')
+
 
 
 def test_CVF2pandas(h2writer):
@@ -73,15 +75,7 @@ def test_CVF2pandas(h2writer):
                  'wind.tint': [0.1, 0.15, 0.2]}
     functions = {'Name': lambda x: 'sim_wsp' + str(x['wind.wsp']) + '_ti' + str(x['wind.tint'])}
 
-    df = h2writer.from_CVF(constants, variables, functions)
-    assert len(df) == 9
-    assert set(list(df)) == set(['simulation.time_stop', 'wind.wsp', 'wind.tint', 'Name'])
-
-
-def test_from_definition(h2writer):
-    from wetb.hawc2.tests.test_files import dlc_definition
-
-    definition_file = dlc_definition.__file__
-    df = h2writer.from_definition(definition_file)
-    assert len(df) == 9
-    assert set(list(df)) == set(['simulation.time_stop', 'wind.wsp', 'wind.tint', 'Name'])
+    h2writer.from_CVF(constants, variables, functions)
+    h2writer.write_all(path + 'htc/tmp/')
+    assert len(h2writer.contents) == 9
+    assert set(list(h2writer.contents)) == set(['simulation.time_stop', 'wind.wsp', 'wind.tint', 'Name'])
diff --git a/wetb/hawc2flow/Evaluate.py b/wetb/hawc2flow/Evaluate.py
new file mode 100644
index 0000000000000000000000000000000000000000..38798fc3b19066a2b1312e1c2797bcf9dba0b88e
--- /dev/null
+++ b/wetb/hawc2flow/Evaluate.py
@@ -0,0 +1,360 @@
+import pandas as pd
+import numpy as np
+import re, os
+import click
+from scipy import signal
+from functools import wraps 
+from wetb.hawc2.Hawc2io import ReadHawc2
+from wetb.fatigue_tools.fatigue import eq_load
+
+    
+def isFloat(string):
+    try:
+        float(string)
+        return True
+    except ValueError:
+        return False
+
+
+class HAWC2DataFrame(pd.DataFrame):
+    """
+    A subclass of a pandas.DataFrame. The HAWC2DataFrame is able to link to a
+    directory of HAWC2 result files. The HAWC2DataFrame contains an Accessor,
+    named 'wetb', which contains functions to  populate the HAWC2DataFrame with
+    post-processed statistics, and to aggregate the statistics.
+    """
+
+    _metadata = ['_fields', 'channels', '_directory', '_filenames']
+    
+    @property
+    def _constructor(self):
+        return HAWC2DataFrame
+        
+        
+    def __init__(self, *args, dir=None, pattern=None, channels=None, **kwargs):
+        
+        if all(x is not None for x in [dir, channels]):            
+            if pattern is None:
+                res = self.link_directory(dir, channels)
+            else:
+                res = self.link_results(dir, pattern, channels)
+            super(HAWC2DataFrame, self).__init__(res)
+
+        else:
+            super(HAWC2DataFrame, self).__init__(*args, **kwargs)
+
+
+    
+    def __finalize__(self, other, method=None, **kwargs):
+        """propagate metadata from other to self """
+        # merge operation: using metadata of the left object
+        if method == 'merge':
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other.left, name, None))
+        # concat operation: using metadata of the first object
+        elif method == 'concat':
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other.objs[0], name, None))
+        else:
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other, name, None))
+        return self
+
+    
+    def link_results(self, directory, pattern_string, channels):
+        self._directory = directory
+        pattern, self._fields = self._compile_pattern(pattern_string)
+
+        # get all filenames that fit the pattern
+        self._filenames = [x[:-4] for x in os.listdir(directory) if x.endswith('.sel')]
+        self._filenames = [x for x in self._filenames if pattern.match(x)]
+        
+        self.channels = channels
+        # Extract input attributes and put in dataframe
+        dat = []
+        for fn in self._filenames:
+            dat.append([float(x) if isFloat(x) else x for x in pattern.findall(fn)[0]])
+        dat = pd.DataFrame(dat)
+        
+        # set column multi index
+        if not dat.empty:
+            column_tuples = list(zip(*[self._fields, ['']*len(self._fields)]))
+            dat.columns = pd.MultiIndex.from_tuples(column_tuples, names=['channel', 'stat'])
+        return dat
+    
+    
+    def link_directory(self, directory, channels):
+        '''
+        Links all result files in a directory. Used if no pattern is given.
+        '''
+        self._directory = directory
+        self.channels = channels
+
+        self._fields = ['filename']
+        # get all filenames of result files
+        self._filenames = [x[:-4] for x in os.listdir(directory) if x.endswith('.sel')]
+        dat = pd.DataFrame(self._filenames)
+
+        # set column multi index
+        if not dat.empty:
+            column_tuples = list(zip(*[self._fields, ['']*len(self._fields)]))
+            dat.columns = pd.MultiIndex.from_tuples(column_tuples, names=['channel', 'stat'])
+        return dat
+        
+            
+    @staticmethod
+    def _compile_pattern(pattern_string):
+        brackets = re.compile('{(.*?)}')
+        fields = brackets.findall(pattern_string)
+        for field in fields:
+            pattern_string = pattern_string.replace('{'+ field +'}', '(.*)')
+        return re.compile(pattern_string), fields   
+        
+                 
+    def __call__(self, **kwargs):
+        return self[self._mask( **kwargs)]
+
+
+    def _mask(self, **kwargs):
+        '''
+        Returns a mask for refering to a dataframe, or self.Data, or self.Data_f, etc.
+        example. dlc.mask(wsp=[12, 14], controller='noIPC')
+        '''
+        N = len(self)
+        mask = [True] * N
+        for key, value in kwargs.items():
+            if isinstance(value, (list, tuple, np.ndarray)):
+                mask_temp = [False] * N
+                for v in value:
+                    mask_temp = mask_temp | (self[key] == v)
+                mask = mask & mask_temp
+            else: #scalar, or single value
+                mask = mask & (self[key] == value)
+        return mask
+
+
+
+@pd.api.extensions.register_dataframe_accessor("wetb")
+class wetbAccessor(object):
+    def __init__(self, pandas_obj):
+        self._validate(pandas_obj)
+        self._obj = pandas_obj
+    
+        
+    @staticmethod
+    def _validate(obj):
+        # verify that a HAWC2DataFrame is using the wetb accessor.
+        if not isinstance(obj, HAWC2DataFrame):
+            raise TypeError(f"wetb namespace is not accessible by type {type(obj)}. Expected: { type(HAWC2DataFrame())}")
+            
+            
+    @classmethod    
+    def populate_method(cls, method_name, label=None):
+        label = label or method_name
+        def decorator(func):
+            @wraps(func) 
+            def wrapper(self, channels, *args, **kwargs): 
+                return self._add_stat(func, label, channels, *args, **kwargs)
+            setattr(cls, method_name, wrapper)
+            # Note we are not binding func, but wrapper which accepts self but does exactly the same as func
+            return func # returning func means func can still be used normally
+        return decorator   
+        
+                 
+    def fetch_result(self, idx, channels=None):
+        '''
+        Returns a time series dataframe of a single result file given an index number.
+        '''
+        fn = self._obj._filenames[idx]
+        channels = channels or self.channels
+
+        raw = ReadHawc2(os.path.join(self._obj._directory, fn)).ReadAll(ChVec=[i-1 for i in channels.values()])
+        #raw = readHawc2Res(os.path.join(self._directory, fn), channels)
+        return pd.DataFrame(raw, columns=channels.keys())
+        # except:
+        #     print(f'File {fn} could not be loaded.')
+        
+            
+    def iter_sim(self, **kwargs):
+        '''
+        Iterates over simulation result files given the input filter defined in kwargs.
+        '''
+        sims_to_iterate = self._obj(**kwargs)
+        for idx, row in sims_to_iterate.iterrows():
+            raw = self.fetch_result(idx)
+            yield row[self._obj._fields], raw
+
+    
+    def _add_stat(self, func, stat_name, channels=None, *args, **kwargs):
+        '''
+        Adds a column of statistics for the given channels using the given function.
+        The function should take a pandas series (1d array) and return a float.
+        '''
+        values = []
+        if channels is None:
+            channels = self._obj.channels
+        else:
+            channels = {k:v for k,v in self._obj.channels.items() if k in channels}
+        
+        channel_string = ', '.join(channels)
+        print(f'Calculating {stat_name} for {channel_string}...')
+        
+        with click.progressbar(self._obj.index) as bar:
+            for idx in bar:
+                raw = self.fetch_result(idx, channels)
+                values.append(func(raw, *args, **kwargs))
+
+        df = pd.DataFrame(values)
+        # add multi index columns
+        col_ch     = list(channels.keys())
+        col_stat   = [stat_name]*len(col_ch)
+        col_tuples = list(zip(*[col_ch, col_stat]))
+        df.columns = pd.MultiIndex.from_tuples(col_tuples, names=['channel', 'stat'])
+
+        return self._obj.join(df)
+
+
+    def aggregate_rows(self, key):
+        # used for taking the mean over all seeds
+        print(f'Calculating mean over key={key}...')
+        in_fields = [x for x in self._obj._fields if x != key]
+        out_fields = [x for x in list(self._obj) if x[0] not in self._obj._fields]
+        in_atts = self._obj[in_fields].drop_duplicates()
+        new_dat = []
+        with click.progressbar(in_atts.iterrows()) as bar:
+            for _, x in bar:
+                filt = {k[0]: v for k, v in dict(x).items()}
+                new_dat.append(list(x) + list(self._obj(**filt)[out_fields].mean().values))
+            
+        # Make new HAWC2DataFrame and copy metadata (is there a better way?)
+        new_df = HAWC2DataFrame(new_dat)
+        for attr in self._obj._metadata:
+            setattr(new_df, attr, getattr(self._obj, attr))
+        new_df._fields = in_fields
+
+        # add multi index columns
+        col_ch     = in_fields + [x[0] for x in out_fields]
+        col_stat   = ['']*len(in_fields) + [x[1] for x in out_fields]
+        col_tuples = list(zip(*[col_ch, col_stat]))
+        new_df.columns = pd.MultiIndex.from_tuples(col_tuples,
+                    names=['channel', 'stat'])
+        return new_df
+
+
+    def aggregate_columns(self, source, dest, drop=False):
+        # used for taking the mean over all blades.
+        channel_string = ', '.join(source)
+        print(f'Calculating mean over {channel_string}...')
+        new_df = self._obj.copy()
+        # get all unique stat indices
+        col_stat = list(set(x[1] for x in list(self._obj)))
+        for stat in col_stat:
+            keys = [x for x in list(self._obj) if x[0] in source and x[1] == stat]
+            if not keys:
+                continue
+            new_df[(dest, stat)] = self._obj[keys].mean(axis=1)
+            
+        if drop:
+            # delete source columns
+            new_df = new_df.drop(columns=source)
+        return new_df
+
+    
+    @staticmethod
+    def read_csv(filename):
+        '''
+        Reads a postproc csv file which was generated by the HAWC2Res class. Returns
+        a DataFrame
+        '''
+        df           = HAWC2DataFrame(pd.read_csv(filename, header =[0, 1]))
+        multicol     = list(df)
+        multicol_new = []
+        for col in multicol:
+            if 'Unnamed' in col[1]:
+                multicol_new.append((col[0], ''))
+            else:
+                multicol_new.append(col)
+        df.columns = pd.MultiIndex.from_tuples(multicol_new)
+        return df
+
+
+@wetbAccessor.populate_method('DEL')
+def DEL(x, m=4):
+    '''
+    Calculates the damage equivalent load of a set of time series loads.
+    args:
+        x (pd.DataFrame): load time series, where each column is a different load channel
+        m (float): wohler constant (default: 4)
+    returns:
+        DEL (list of floats): Damage equivalent load for each load channel.
+    '''
+    DEL = []
+    for k in list(x):
+        DEL.append(eq_load(x[k].values, m=m)[0][0])
+    return DEL
+    
+    
+@wetbAccessor.populate_method('mean', label='Mean')
+def _mean(x):
+    '''
+    Calculates the mean of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): Mean of each time series.
+    '''
+    return x.mean().values
+    
+    
+@wetbAccessor.populate_method('var', label='Var')
+def _var(x):
+    '''
+    Calculates the variance of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): variance of each time series.
+    '''
+    return x.var().values
+
+
+@wetbAccessor.populate_method('std', label='Std')
+def _std(x):
+    '''
+    Calculates the standard deviation of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): standard deviation of each time series.
+    '''
+    return x.std().values
+
+
+@wetbAccessor.populate_method('final')
+def _final(x):
+    '''
+    Returns the final value of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): final value of each time series.
+    '''
+    return x.iloc[-1].values
+
+
+@wetbAccessor.populate_method('PSD')
+def _psd(x, fs=100):
+    '''
+    Calculates power spectral density (PSD) of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+        fs (float): sampling frequency in Hz (default: 100Hz).
+    returns:
+        (list of arrays): PSD of each time series.
+    '''
+    PSD = []
+    for k in list(x):
+        f, Py = signal.welch(x[k].values, fs=fs, nperseg=1024*8)
+        PSD.append(Py)
+    return PSD
+            
diff --git a/wetb/hawc2flow/Generate.py b/wetb/hawc2flow/Generate.py
new file mode 100644
index 0000000000000000000000000000000000000000..7137b7abb886900eaf48496201093bd7c0dc903e
--- /dev/null
+++ b/wetb/hawc2flow/Generate.py
@@ -0,0 +1 @@
+from wetb.hawc2.hawc2_input_writer import HAWC2InputWriter, JinjaWriter
diff --git a/wetb/hawc2flow/Simulate.py b/wetb/hawc2flow/Simulate.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/hawc2flow/__init__.py b/wetb/hawc2flow/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391