Repository: NREL/bifacial_radiance (language: python)

Issue:
Is it possible to pass custom sensors like tutorial 12 in a 1-axis analysis?
Hi again Bifacial Radiance team!

As my post title hints, I'm trying to make ground irradiance measurements like the ones in the APV tutorial 12 but using the 1-axis tracking workflow. As far as I can tell, the scan locations for the front and back of the PV surfaces are passed through the trackerdict, so the moduleAnalysis function isn't applicable here, or? Analysis1axis doesn't seem to be capable of this out of the box either.

Is there a way to change scan locations in the trackerdict? Is there potentially a way to append a custom scan location to 1axis scenes? If my questions are waaay far off, I hope you understand the problem statement and I'd be grateful for any help :)

Hope you've had a nice weekend.
Back and front Analysis separation

gendaylit1axis not finding match on EPW
EPW files are being shifted 1 hour to match TMY data. However this creates 1 extra entry the next year.
![image](https://user-images.githubusercontent.com/22824806/106505740-f47a0f00-6485-11eb-90cd-995ecfdc75c6.png)



![image](https://user-images.githubusercontent.com/22824806/106505719-ecba6a80-6485-11eb-879e-f70b544d6136.png)


The current trackerdict matching between wanted startdate and enddate is very archaic, and it's failing to find 01/01 correctly since it ifnds it in the next year. 

![image](https://user-images.githubusercontent.com/22824806/106505796-0956a280-6486-11eb-9d4b-e3c9373e1c8f.png)

![image](https://user-images.githubusercontent.com/22824806/106505882-29866180-6486-11eb-98b7-aff3336617ca.png)



improve addMaterials to Ground.rad
addMaterial to Ground.rad function currently just adds new materials at the end, and does not double-check if material exists and if it needs to be overwritten, so it ends up being unclear which material if they are repeated names gets used. 


gendaylit1axis not finding match on EPW
EPW files are being shifted 1 hour to match TMY data. However this creates 1 extra entry the next year.
![image](https://user-images.githubusercontent.com/22824806/106505740-f47a0f00-6485-11eb-90cd-995ecfdc75c6.png)



![image](https://user-images.githubusercontent.com/22824806/106505719-ecba6a80-6485-11eb-879e-f70b544d6136.png)


The current trackerdict matching between wanted startdate and enddate is very archaic, and it's failing to find 01/01 correctly since it ifnds it in the next year. 

![image](https://user-images.githubusercontent.com/22824806/106505796-0956a280-6486-11eb-9d4b-e3c9373e1c8f.png)

![image](https://user-images.githubusercontent.com/22824806/106505882-29866180-6486-11eb-98b7-aff3336617ca.png)



GenCumSky cannot read metdata weather
In development branch gencumsky() fails to properly read tmy3_temp.csv or epw_temp.csv.
using `bifacial_radiance.readWeatherFile()` on the following files:

```
weatherFile = 'bifacial_radiance\tests\USA_CO_Boulder.724699_TMY2.epw'
radObj.readWeatherFile(weatherFile,label='center')
radObj.gencumsky()
```
![image](https://user-images.githubusercontent.com/75647238/121060520-00514700-c780-11eb-97ba-772775e09b44.png)

```
weatherFile = 'bifacial_radiance\tests\724666TYA.CSV'
radObj.readWeatherFile(weatherFile,label='center')
radObj.gencumsky()
```
![image](https://user-images.githubusercontent.com/75647238/121060617-1fe86f80-c780-11eb-9a64-de205aa2e845.png)


GenCumSky can successfully read EPW or CSV directly as below
```
weatherFile = 'bifacial_radiance\tests\USA_CO_Boulder.724699_TMY2.epw'
radObj.gencumsky(weatherFile)
```
![image](https://user-images.githubusercontent.com/75647238/121062941-f2e98c00-c782-11eb-8ac3-b97cba28aa1c.png)


```
weatherCSV = 'bifacial_radiance\tests\724666TYA.CSV'
radObj.gencumsky(weatherCSV)
```
![image](https://user-images.githubusercontent.com/75647238/121063077-1ca2b300-c783-11eb-8240-c49338a0cb4e.png)



Discussion / hints:
Hello :)

Yes, that makes total sense, and yes, we need to add that. I think I just pushed some changes into the a development-version subbranch `modScan-Analysis1axis` that should do what you want, with the caveat that:

1) You must pass a frontscan and/or backscan dictionary that have the right parameters 
2) Frontscan and backscan (if you pass both, or if you pass 1 with the generated) must have the same number of sensors, else it will complain when saving results.

This mod will only sample the same spot for all times simulated though, will not vary with the angle, so it's ideal for sampling the ground for eample for the APV. Do you want to still scan on the module but shifted up and down? If that's the case on the lines I modified we have to pass ability to modify parameters of the calculated frontscan/backscan, instead of a whole new frontscan/backscan... Let me know. We have to think how to do both elegantly without adding too many parameters.
 
FYI The code I added is in ``analysis1axis`` function:

                frontscan, backscan = analysis.moduleAnalysis(scene=scene, modWanted=modWanted, rowWanted=rowWanted, sensorsy=sensorsy)
                if modscanfront is not None:
                         frontscan = modscanfront
                if modscanback is not None:
                        backscan = modscanback



@metavektor my last push into the branch allows you to change any/all of the keys in the frontscan backscan to analysis1axis. Hope this helps! Closing the issue, this feature will be officially released in next release


Hey hey, wow you're too fast for me! :) 

I ended up fudging a solution for my weird problem and essentially scripted static scenes with tilts read in through a for loop, which is definitely less than ideal but physically the same. It definitely sounds like you've added a functionality that will help avoid this! On the weekend I'll jump back in to test your addition out.

Thank you again for your time and awesome efforts!






Fix patch:
diff --git a/MANIFEST.in b/MANIFEST.in
index 2471e239..ca7e2197 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,2 +1,4 @@
 include versioneer.py
 include bifacial_radiance/_version.py
+include bifacial_radiance/data/*
+include bifacial_radiance/images/*
\ No newline at end of file
diff --git a/README.md b/README.md
index 32030478..a64911d8 100644
--- a/README.md
+++ b/README.md
@@ -70,7 +70,6 @@ If you use bifacial_radiance in a published work, please cite:
 
     Ayala Pelaez, Silvana, Deline, Chris (2020). bifacial_radiance: a python package for modeling bifacial solar photovoltaic systems. Journal of Open Source Software, 5(50), 1865, https://doi.org/10.21105/joss.01865
     
-
 Additional bifacial_radiance publications with validation of the software include:
 * Deline, Chris, and Ayala, Silvana. Bifacial_Radiance. Computer Software. https://github.com/NREL/bifacial_radiance. USDOE Office of     Energy Efficiency and Renewable Energy (EERE), Solar Energy Technologies Office (EE-4S). 17 Dec. 2017. Web. doi:10.11578/dc.20180530.16. https://www.osti.gov/doecode/biblio/6869
 * Ayala Pelaez S, Deline C, Greenberg P, Stein JS, Kostuk RK. Model and validation of single-axis tracking with bifacial PV. IEEE J Photovoltaics. 2019;9(3):715-721. https://ieeexplore.ieee.org/document/8644027 and https://www.nrel.gov/docs/fy19osti/72039.pdf (pre-print, conference version)
diff --git a/bifacial_radiance/__init__.py b/bifacial_radiance/__init__.py
index 27e99daa..55146867 100644
--- a/bifacial_radiance/__init__.py
+++ b/bifacial_radiance/__init__.py
@@ -4,6 +4,7 @@
 from bifacial_radiance import modelchain
 from bifacial_radiance.gui import gui
 from bifacial_radiance import mismatch
+from bifacial_radiance import spectral_utils
 from ._version import get_versions
 __version__ = get_versions()['version']
 del get_versions
diff --git a/bifacial_radiance/data/module.json b/bifacial_radiance/data/module.json
index 3ace8d96..a830f799 100644
--- a/bifacial_radiance/data/module.json
+++ b/bifacial_radiance/data/module.json
@@ -77,22 +77,24 @@
         "bifi": 0.9,
         "cellModule": null,
         "modulefile": "objects\\Prism_Solar_Bi60.rad",
-        "moduleoffset": 0.15,
-        "numpanels": 1,
-        "scenex": 1.084,
-        "sceney": 1.695,
+        "modulematerial": "black",
+        "numpanels": 2,
+        "offsetfromaxis": 0.15,
+        "scenex": 1.03,
+        "sceney": 4.11,
         "scenez": 0.15,
-        "text": "! genbox black Prism_Solar_Bi60 0.984 1.695 0.02 | xform -t -0.492 -0.8475 0.15 -a 1 -t 0 1.695 0",
+        "text": "! genbox black Prism_Solar_Bi60 0.98 1.98 0.02 | xform -t -0.49 -2.055 0.15 -a 2 -t 0 2.13 0\r\n! genrev Metal_Grey tube1 t*1.03 0.05 32 | xform -ry 90 -t -0.515 0 0",
         "torquetube": {
-            "bool": false,
+            "bool": true,
             "diameter": 0.1,
             "material": "Metal_Grey",
             "tubetype": "round"
         },
-        "x": 0.984,
-        "xgap": 0.1,
-        "y": 1.695,
-        "ygap": 0.0,
+        "x": 0.98,
+        "xgap": 0.05,
+        "y": 1.98,
+        "ygap": 0.15,
+        "z": 0.02,
         "zgap": 0.1
     },
     "Prism Solar Bi60 landscape": {
diff --git a/bifacial_radiance/main.py b/bifacial_radiance/main.py
index bef51a61..ccb32de3 100644
--- a/bifacial_radiance/main.py
+++ b/bifacial_radiance/main.py
@@ -122,6 +122,30 @@ def _interactive_directory(title=None):
     root.attributes("-topmost", True) #Bring to front
     return filedialog.askdirectory(parent=root, title=title)
 
+def _modDict(originaldict, moddict):
+    '''
+    Compares keys in originaldict with moddict and updates values of 
+    originaldict to moddict if existing.
+    
+    Parameters
+    ----------
+    originaldict : dictionary
+        Original dictionary calculated, for example frontscan or backscan dictionaries.
+    moddict : dictionary
+        Modified dictionary, for example modscan['xstart'] = 0 to change position of x.
+    
+    Returns
+    -------
+    originaldict : dictionary
+        Updated original dictionary with values from moddict.
+    '''
+    for key in moddict:
+        try:
+            originaldict[key] = moddict[key]
+        except:
+            print("Wrong key in modified dictionary")
+                
+    return originaldict
 
 class RadianceObj:
     """
@@ -262,25 +286,51 @@ def save(self, savefile=None):
             pickle.dump(self, f)
         print('Saved to file {}'.format(savefile))
 
-    def addMaterial(self, material, Rrefl, Grefl, Brefl, materialtype='plastic', spec=0, rough=0, material_file=None, comment=None):
+    def addMaterial(self, material, Rrefl, Grefl, Brefl, materialtype='plastic', spec=0, rough=0, material_file=None, comment=None, rewrite=True):
     
         if material_file is None:
             material_file = 'ground.rad'    
     
         matfile = os.path.join('materials', material_file)
-    
-        # append -- This will create the file if it doesn't exist
-        file_object = open(matfile, 'a')  
-        file_object.write("\n\n")
-        if comment is not None:
-            file_object.write("#{}".format(comment))
-        file_object.write("\nvoid {} {}".format(materialtype, material))
-        if materialtype == 'glass':
-            file_object.write("\n0\n0\n3 {} {} {}".format(Rrefl, Grefl, Brefl))
-        else:
-            file_object.write("\n0\n0\n5 {} {} {} {} {}".format(Rrefl, Grefl, Brefl, spec, rough))
-        file_object.close()
-        print('Added material {} to file {}'.format(material, material_file))
+        
+        with open(matfile, 'r') as fp:
+            buffer = fp.readlines()
+                
+        # search buffer for material matching requested addition
+        found = False
+        for i in buffer:
+            if materialtype and material in i:
+                loc = buffer.index(i)
+                found = True
+                break
+        if found:
+            if rewrite:            
+                print('Material exists, overwriting...\n')
+                if comment is None:
+                    pre = loc - 1
+                else:
+                    pre = loc - 2            
+                # commit buffer without material match
+                with open(matfile, 'w') as fp:
+                    for i in buffer[0:pre]:
+                        fp.write(i)
+                    for i in buffer[loc+4:]:
+                        fp.write(i)
+        if (found and rewrite) or (not found):
+            # append -- This will create the file if it doesn't exist
+            file_object = open(matfile, 'a')
+            file_object.write("\n\n")
+            if comment is not None:
+                file_object.write("#{}".format(comment))
+            file_object.write("\nvoid {} {}".format(materialtype, material))
+            if materialtype == 'glass':
+                file_object.write("\n0\n0\n3 {} {} {}".format(Rrefl, Grefl, Brefl))
+            else:
+                file_object.write("\n0\n0\n5 {} {} {} {} {}".format(Rrefl, Grefl, Brefl, spec, rough))
+            file_object.close()
+            print('Added material {} to file {}'.format(material, material_file))
+        if (found and not rewrite):
+            print('Material already exists\n')
 
     def exportTrackerDict(self, trackerdict=None,
                           savefile=None, reindex=None):
@@ -700,7 +750,7 @@ def _convertTMYdate(data, meta):
             endtime = dd[0]+'_'+dd[1] + '_23'
         
         tmydata_trunc = self._saveTempTMY(tmydata,'tmy3_temp.csv', 
-                                              starttime=starttime, endtime=endtime)
+                                          starttime=starttime, endtime=endtime)
 
         if daydate is not None:  # also remove GHI = 0 for HPC daydate call.
             tmydata_trunc = tmydata_trunc[tmydata_trunc.GHI > 0]
@@ -948,7 +998,6 @@ def gendaylit(self, timeindex, metdata=None, debug=False):
         time = metdata.datetime[timeindex]
         #filename = str(time)[2:-9].replace('-','_').replace(' ','_').replace(':','_')
         filename = time.strftime('%Y_%m_%d_%H_%M')
-        print(filename)
         skyname = os.path.join(sky_path,"sky2_%s_%s_%s.rad" %(lat, lon, filename))
 
         skyFile = open(skyname, 'w')
@@ -1638,7 +1687,7 @@ def makeModule(self, name=None, x=None, y=None, z=None, bifi=1, modulefile=None,
         # Adding the option to replace the module thickess
         if z is None:
             z = 0.020
-        
+            
         if modulematerial is None:
             modulematerial = 'black'
             
@@ -1655,9 +1704,9 @@ def makeModule(self, name=None, x=None, y=None, z=None, bifi=1, modulefile=None,
                     packagingfactor = 100.0
 
                 except NameError as err: # probably because no x or y passed
-                    raise Exception('makeModule variable {} and cellLevelModule'+
-                                    'Params is None.  One or the other must'+
-                                    ' be specified.'.format(err.args[0]))
+                    raise Exception('makeModule variable {}'.format(err.args[0])+
+                                    ' and cellLevelModuleParams is None.  '+
+                                    'One or the other must be specified.')
             else:
                 c = cellLevelModuleParams
                 x = c['numcellsx']*c['xcell'] + (c['numcellsx']-1)*c['xcellgap']
@@ -2237,7 +2286,8 @@ def makeScene1axis(self, trackerdict=None, moduletype=None, sceneDict=None,
 
 
     def analysis1axis(self, trackerdict=None, singleindex=None, accuracy='low',
-                      customname=None, modWanted=None, rowWanted=None, sensorsy=9, hpc=False):
+                      customname=None, modWanted=None, rowWanted=None, sensorsy=9, hpc=False,
+                      modscanfront = None, modscanback = None):
         """
         Loop through trackerdict and runs linescans for each scene and scan in there.
 
@@ -2257,6 +2307,22 @@ def analysis1axis(self, trackerdict=None, singleindex=None, accuracy='low',
             Row to be sampled. Index starts at 1. (row 1)
         sensorsy : int 
             Sampling resolution for the irradiance across the collector width.
+        modscanfront : dict
+            dictionary with one or more of the following key: xstart, ystart, zstart, 
+            xinc, yinc, zinc, Nx, Ny, Nz, orient. All of these keys are ints or 
+            floats except for 'orient' which takes x y z values as string 'x y z'
+            for example '0 0 -1'. These values will overwrite the internally
+            calculated frontscan dictionary for the module & row selected. If modifying 
+            Nx, Ny or Nz, make sure to modify on modscanback to avoid issues on 
+            results writing stage. 
+        modscanback : dict
+            dictionary with one or more of the following key: xstart, ystart, zstart, 
+            xinc, yinc, zinc, Nx, Ny, Nz, orient. All of these keys are ints or 
+            floats except for 'orient' which takes x y z values as string 'x y z'
+            for example '0 0 -1'. These values will overwrite the internally
+            calculated frontscan dictionary for the module & row selected.  If modifying 
+            Nx, Ny or Nz, make sure to modify on modscanback to avoid issues on 
+            results writing stage. 
 
         Returns
         -------
@@ -2305,7 +2371,9 @@ def analysis1axis(self, trackerdict=None, singleindex=None, accuracy='low',
             try:  # look for missing data
                 analysis = AnalysisObj(octfile,name)
                 name = '1axis_%s%s'%(index,customname,)
-                frontscan, backscan = analysis.moduleAnalysis(scene=scene, modWanted=modWanted, rowWanted=rowWanted, sensorsy=sensorsy)
+                frontscan, backscan = analysis.moduleAnalysis(scene=scene, modWanted=modWanted, 
+                                                rowWanted=rowWanted, sensorsy=sensorsy, 
+                                                modscanfront=modscanfront, modscanback=modscanback)
                 analysis.analysis(octfile=octfile,name=name,frontscan=frontscan,backscan=backscan,accuracy=accuracy, hpc=hpc)                
                 trackerdict[index]['AnalysisObj'] = analysis
             except Exception as e: # problem with file. TODO: only catch specific error types here.
@@ -2352,7 +2420,8 @@ def analysis1axis(self, trackerdict=None, singleindex=None, accuracy='low',
             if self.cumulativesky is True: 
                 frontcum = pd.DataFrame()
                 rearcum = pd.DataFrame()
-                temptrackerdict = trackerdict[0.0]['AnalysisObj']
+                temptrackerdict = trackerdict[list(trackerdict)[0]]['AnalysisObj']
+                #temptrackerdict = trackerdict[0.0]['AnalysisObj']
                 frontcum ['x'] = temptrackerdict.x
                 frontcum ['y'] = temptrackerdict.y
                 frontcum ['z'] = temptrackerdict.z
@@ -2374,7 +2443,9 @@ def analysis1axis(self, trackerdict=None, singleindex=None, accuracy='low',
                     cumscene.sceneDict['tilt']=0
                     cumscene.sceneDict['clearance_height'] = self.hub_height
                     cumanalysisobj = AnalysisObj()
-                    frontscan, backscan = cumanalysisobj.moduleAnalysis(scene=cumscene, modWanted=modWanted, rowWanted=rowWanted, sensorsy = sensorsy)
+                    frontscan, backscan = cumanalysisobj.moduleAnalysis(scene=cumscene, modWanted=modWanted, 
+                                                rowWanted=rowWanted, sensorsy=sensorsy, 
+                                                modscanfront=modscanfront, modscanback=modscanback)
                     x,y,z = cumanalysisobj._linePtsArray(frontscan)
                     x,y,rearz = cumanalysisobj._linePtsArray(backscan)
         
@@ -2448,7 +2519,7 @@ def __init__(self, materialOrAlbedo=None, material_file=None):
             
         self.material_file = material_file           
         if materialOrAlbedo is None: # Case where it's none.
-            print('\nInput albedo 0-1, or string from GroundObj.printGroundMaterials().'
+            print('\nInput albedo 0-1, or string from ground.printGroundMaterials().'
             '\nAlternatively, run setGround after readWeatherData()'
             'and setGround will read metdata.albedo if availalbe')
             return
@@ -2535,8 +2606,8 @@ def printGroundMaterials(self, materialString=None):
             try:
                 index = _findme(keys,materialString)[0]
             except IndexError:
-                warnings.warn(f'Error - materialString not in '
-                              '{self.material_file}: {materialString}')
+                warnings.warn('Error - materialString not in '
+                              f'{self.material_file}: {materialString}')
             return(np.array([[Rreflall[index], Greflall[index], Breflall[index]]]))
         else:
             return(keys)
@@ -3007,6 +3078,17 @@ def __init__(self, tmydata, metadata, label = 'right'):
         self.dni = np.array(tmydata.DNI)
         self.albedo = np.array(tmydata.Alb)
         
+        # Try and retrieve dewpoint and pressure
+        try:
+            self.dewpoint = np.array(tmydata['temp_dew'])
+        except KeyError:
+            self.dewpoint = None
+
+        try:
+            self.pressure = np.array(tmydata['atmospheric_pressure'])
+        except KeyError:
+            self.pressure = None
+
         #v0.2.5: initialize MetObj with solpos, sunrise/set and corrected time
         datetimetz = pd.DatetimeIndex(self.datetime)
         try:  # make sure the data is tz-localized.
@@ -3608,9 +3690,9 @@ def _saveResults(self, data, reardata=None, savefile=None, RGB = False):
         
         if RGB:
             data_sub = {key:data[key] for key in ['x', 'y', 'z', 'r', 'g', 'b', 'Wm2', 'mattype']}
-            self.R = data['R']
-            self.G = data['G']
-            self.B = data['B']
+            self.R = data['r']
+            self.G = data['g']
+            self.B = data['b']
             self.x = data['x']
             self.y = data['y']
             self.z = data['z']
@@ -3645,9 +3727,10 @@ def _saveResults(self, data, reardata=None, savefile=None, RGB = False):
                 data_sub['rearB'] = self.rearB
                 
                 df = pd.DataFrame.from_dict(data_sub)
+                df.reindex(columns=['x','y','z','rearZ','mattype','rearMat',
+                                    'Wm2Front','Wm2Back','Back/FrontRatio',
+                                    'R','G','B', 'rearR','rearG','rearB'])
                 df.to_csv(os.path.join("results", savefile), sep = ',',
-                          columns = ['x','y','z','rearZ','mattype','rearMat',
-                                     'Wm2Front','Wm2Back','Back/FrontRatio', 'R','G','B', 'rearR','rearG','rearB'],
                                      index = False) # new in 0.2.3
 
             else:
@@ -3719,7 +3802,8 @@ def _saveResultsCumulative(self, data, reardata=None, savefile=None):
         return (savefile)
 
     def moduleAnalysis(self, scene, modWanted=None, rowWanted=None,
-                       sensorsy=9.0, frontsurfaceoffset=0.001, backsurfaceoffset=0.001, debug=False):
+                       sensorsy=9.0, frontsurfaceoffset=0.001, backsurfaceoffset=0.001, 
+                       modscanfront=None, modscanback=None, debug=False):
         """
         This function defines the scan points to be used in the 
         :py:class:`~bifacial_radiance.AnalysisObj.analysis` function,
@@ -3738,7 +3822,25 @@ def moduleAnalysis(self, scene, modWanted=None, rowWanted=None,
             (CW) of the module(s)
         debug : bool
             Activates various print statemetns for debugging this function.
-
+        modscanfront : dict
+            Dictionary to modify the fronstcan values established by this routine 
+            and set a specific value. Keys possible are 'xstart', 'ystart', 'zstart',
+            'xinc', 'yinc', 'zinc', 'Nx', 'Ny', 'Nz', and 'orient'. If modifying 
+            Nx, Ny or Nz, make sure to modify on modscanback to avoid issues on 
+            results writing stage. All of these keys are ints or 
+            floats except for 'orient' which takes x y z values as string 'x y z'
+            for example '0 0 -1'. These values will overwrite the internally
+            calculated frontscan dictionary for the module & row selected.
+        modscanback: dict
+            Dictionary to modify the backscan values established by this routine 
+            and set a specific value. Keys possible are 'xstart', 'ystart', 'zstart',
+            'xinc', 'yinc', 'zinc', 'Nx', 'Ny', 'Nz', and 'orient'. If modifying 
+            Nx, Ny or Nz, make sure to modify on modscanback to avoid issues on 
+            results writing stage. All of these keys are ints or 
+            floats except for 'orient' which takes x y z values as string 'x y z'
+            for example '0 0 -1'. These values will overwrite the internally
+            calculated frontscan dictionary for the module & row selected.
+        
         Returns
         -------
         frontscan : dictionary
@@ -3943,7 +4045,7 @@ def moduleAnalysis(self, scene, modWanted=None, rowWanted=None,
             print("Offset Shift", x3, y3, z3)
             print("Final Start Coordinate Front", xstartfront, ystartfront, zstartfront)
             print("Increase Coordinates", xinc, yinc, zinc)
-            
+        
         frontscan = {'xstart': firstsensorxstartfront, 'ystart': firstsensorystartfront,
                      'zstart': firstsensorzstartfront,
                      'xinc':xinc, 'yinc': yinc,
@@ -3953,9 +4055,15 @@ def moduleAnalysis(self, scene, modWanted=None, rowWanted=None,
                      'xinc':xinc, 'yinc': yinc,
                      'zinc':zinc, 'Nx': 1, 'Ny':sensorsy, 'Nz':1, 'orient':back_orient }
 
+        if modscanfront is not None:
+            frontscan = _modDict(frontscan, modscanfront)
+        if modscanback is not None:
+            backscan = _modDict(backscan, modscanback)
+                    
         return frontscan, backscan
 
-    def analysis(self, octfile, name, frontscan, backscan, plotflag=False, accuracy='low', hpc = False):
+    def analysis(self, octfile, name, frontscan, backscan,
+                 plotflag=False, accuracy='low', RGB=False, hpc=False):
         """
         General analysis function, where linepts are passed in for calling the
         raytrace routine :py:class:`~bifacial_radiance.AnalysisObj._irrPlot` 
@@ -4002,7 +4110,7 @@ def analysis(self, octfile, name, frontscan, backscan, plotflag=False, accuracy=
                                    plotflag=plotflag, accuracy=accuracy, hpc = hpc)
         # don't save if _irrPlot returns an empty file.
         if frontDict is not None:
-            self._saveResults(frontDict, backDict,'irr_%s.csv'%(name) )
+            self._saveResults(frontDict, backDict,'irr_%s.csv'%(name), RGB=RGB)
 
         return frontDict, backDict
     
@@ -4113,57 +4221,3 @@ def quickExample(testfolder=None):
 
     return analysis
 
-
-if __name__ == "__main__":
-    '''
-    Example of how to run a Radiance routine for a simple rooftop bifacial system
-
-    '''
-    import multiprocessing as mp
-        
-    #  print("This is daydate %s" % (daydate))
-    demo = RadianceObj(simulationname,path=testfolder)
-    #epwfile = demo.getEPW(44, -110)
-    #print(epwfile)
-    demo.setGround(albedo)
-    # metdata = demo.readWeatherFile(epwfile)
-    # moduleDict=demo.makeModule(name=moduletype,x=x,y=y,bifi=bifi, 
-    #                       torquetube=torqueTube, diameter = diameter, tubetype = tubetype, 
-    #                       material = torqueTubeMaterial, zgap = zgap, numpanels = numpanels, ygap = ygap, 
-    #                       rewriteModulefile = True, xgap=xgap, 
-    #                       axisofrotationTorqueTube=axisofrotationTorqueTube)
-    
-    sceneDict = {'module_type':moduletype, 'pitch': pitch, 'hub_height':hub_height, 'nMods':nMods, 'nRows':nRows}  
-    
-    cores = mp.cpu_count()
-    pool = mp.Pool(processes=cores)
-    res = None
-    print ("This is cores", cores)
-
-    try:
-        nodeID = int(os.environ['SLURM_NODEID'])
-    except: 
-        nodeID = 0
-    day_index = (36 * (nodeID))
-    
-    # doing less days for testing
-    start = datetime.datetime.strptime("01-01-2014", "%d-%m-%Y")
-    end = datetime.datetime.strptime("31-12-2014", "%d-%m-%Y") # 2014 not a leap year.
-    #start = datetime.datetime.strptime("14-02-2014", "%d-%m-%Y")
-    #end = datetime.datetime.strptime("26-02-2014", "%d-%m-%Y") # 2014 not a leap year.
-    date_generated = [start + datetime.timedelta(days=x) for x in range(0, (end-start).days)]
-    
-    daylist = []
-    for date in date_generated:
-        daylist.append(date.strftime("%m_%d"))
-    # loop doesn't add last day :
-
-#    print (daylist)
-    for job in range(cores):
-        if day_index+job>=60: #len(daylist):
-            break
-        pool.apply_async(runJob, (daylist[day_index+job],))
-        
-    pool.close()
-    pool.join()
-    pool.terminate()
diff --git a/bifacial_radiance/spectral_utils.py b/bifacial_radiance/spectral_utils.py
new file mode 100644
index 00000000..135b633f
--- /dev/null
+++ b/bifacial_radiance/spectral_utils.py
@@ -0,0 +1,311 @@
+import numpy as np
+import pandas as pd
+from collections.abc import Iterable
+import os
+from scipy import integrate
+
+
+class spectral_property(object):
+    def load_file(filepath):
+        with open(filepath, 'r') as infile:
+            meta = next(infile)[:-1]
+            data = pd.read_csv(infile)
+        
+        return spectral_property(data['value'], data['wavelength'],
+                                 interpolation=meta.split(':')[1])
+    
+    def to_nm(wavelength, units):
+        unit_conversion = { 'nm': 1,
+                            'um': 1000 }
+        
+        # Verify units are in conversion table
+        if units not in unit_conversion:
+            print("Warning: Unknown unit specified. Options are {}.".format(
+                unit_conversion.keys()))
+            units = 'nm'
+            
+        return wavelength * unit_conversion[units]
+    
+    def _linear_interpolation(self, wavelength_nm):
+        # Find upper and lower index
+        upper_bound = self.df[self.df.index > wavelength_nm].index.min()
+        lower_bound = self.df[self.df.index < wavelength_nm].index.max()
+        
+        # Determine values of surrounding indices
+        upper_val = self.df['value'][upper_bound]
+        lower_val = self.df['value'][lower_bound]
+        
+        # Calculate deltas
+        delta_lambda = upper_bound - lower_bound
+        delta_val = upper_val - lower_val
+        
+        return lower_val + delta_val*(wavelength_nm - lower_bound)/delta_lambda
+    
+    def _nearest_interpolation(self, wavelength_nm):
+        # Find upper and lower index
+        upper_bound = self.df[self.df.index > wavelength_nm].index.min()
+        lower_bound = self.df[self.df.index < wavelength_nm].index.max()
+        
+        # Determine which index is closer
+        if (upper_bound - wavelength_nm) < (wavelength_nm - lower_bound):
+            return self.df['value'][upper_bound]
+        return self.df['value'][lower_bound]
+    
+    def _lower_interpolation(self, wavelength_nm):
+        # Find lower index
+        lower_bound = self.df[self.df.index < wavelength_nm].index.max()
+        
+        return self.df['value'][lower_bound]
+    
+    def _upper_interpolation(self, wavelength_nm):
+        # Find upper index
+        upper_bound = self.df[self.df.index > wavelength_nm].index.min()
+        
+        return self.df['value'][upper_bound]
+    
+    interpolation_methods = {
+        'linear':   _linear_interpolation,
+        'nearest':  _nearest_interpolation,
+        'lower':    _lower_interpolation,
+        'upper':    _upper_interpolation
+    }
+    
+    def __init__(self, values, index, index_units='nm', interpolation=None):
+        # Verify lengths match
+        if len(values) != len(index):
+            print("Warning: Length of values and index must match.")
+            return
+        
+        # Convert inputs to list
+        values = [ val for val in values ]
+        index = [ spectral_property.to_nm(idx, index_units) for idx in index ]
+        
+        # Create DataFrame
+        self.df = pd.DataFrame()
+        self.df['value'] = values
+        self.df['wavelength'] = index
+        self.df = self.df.set_index('wavelength')
+        
+        self.interpolation = None
+        if interpolation in spectral_property.interpolation_methods:
+            self.interpolation = \
+                    spectral_property.interpolation_methods[interpolation]
+            self.interpolation_type = interpolation
+        elif interpolation:
+            print("Warning: Specified interpolation type unknown.")
+            
+    def _get_single(self, wavelength, units):
+        # Convert wavelength to nm
+        wavelength = spectral_property.to_nm(wavelength, units)
+        
+        if wavelength in self.df.index:
+            # If the value for that wavelength is known, return it
+            return self.df['value'][wavelength]
+        elif self.interpolation:
+            # Check wavelength is within range
+            if wavelength < self.df.index.min() or \
+               wavelength > self.df.index.max():
+                print("Warning: Requested wavelength outside spectrum.")
+                return None
+            
+            # Return interpolated value
+            return self.interpolation(self, wavelength)
+        
+        return None
+    
+    def __getitem__(self, wavelength, units='nm'):
+        if isinstance(wavelength, Iterable):
+            return np.array([ self._get_single(wl, units) for wl in wavelength ])
+        return self._get_single(wavelength, units)
+    
+    def to_file(self, filepath, append=False):
+        mode = 'w'
+        if append:
+            mode = 'a'
+            
+        with open(filepath, mode) as outfile:
+            outfile.write(f"interpolation:{self.interpolation_type}\n")
+            self.df.to_csv(outfile)
+    
+    def range(self):
+        # Find upper and lower index
+        upper_bound = self.df.index.max()
+        lower_bound = self.df.index.min()
+        
+        return (lower_bound, upper_bound)
+    
+    def scale_values(self, scaling_factor):
+        self.df['value'] *= scaling_factor
+    
+def spectral_albedo_smarts(zen, azm, material, min_wavelength=300,
+                           max_wavelength=4000):
+
+    import pySMARTS
+    
+    smarts_res = pySMARTS.SMARTSSpectraZenAzm('30 31', str(zen), str(azm), material,
+                                     min_wvl=str(min_wavelength),
+                                     max_wvl=str(max_wavelength))
+    
+    return spectral_property(smarts_res['Zonal_ground_reflectance'],
+                             smarts_res['Wvlgth'], interpolation='linear')
+
+def spectral_irradiance_smarts(zen, azm, min_wavelength=300,
+                           max_wavelength=4000):
+
+    import pySMARTS
+
+    smarts_res = pySMARTS.SMARTSSpectraZenAzm('2 3 4', str(zen), str(azm),
+                                     min_wvl=str(min_wavelength),
+                                     max_wvl=str(max_wavelength))
+    
+    dni_spectrum = spectral_property(smarts_res['Direct_normal_irradiance'],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    dhi_spectrum = spectral_property(smarts_res['Difuse_horizn_irradiance'],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    ghi_spectrum = spectral_property(smarts_res['Global_horizn_irradiance'],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    
+    return (dni_spectrum, dhi_spectrum, ghi_spectrum)
+
+
+
+def spectral_irradiance_smarts_SRRL(YEAR, MONTH, DAY, HOUR, ZONE,
+                                LATIT, LONGIT, ALTIT,
+                                RH, TAIR, SEASON, TDAY, SPR, W,
+                                TILT, WAZIM, HEIGHT,
+                                ALPHA1, ALPHA2, OMEGL, GG, BETA, TAU5,
+                                RHOG, material,
+                                IOUT='2 3 4', min_wvl='280', max_wvl='4000'):
+    
+    import pySMARTS
+
+    smarts_res = pySMARTS.SMARTSSRRL(IOUT=IOUT, YEAR=YEAR,MONTH=MONTH,DAY=DAY,HOUR=HOUR, ZONE=ZONE,
+                            LATIT=LATIT, LONGIT=LONGIT, ALTIT=ALTIT, 
+                             RH=RH, TAIR=TAIR, SEASON=SEASON, TDAY=TDAY, SPR=SPR, W=W, 
+                             TILT=TILT, WAZIM=WAZIM, HEIGHT=HEIGHT,
+                             ALPHA1 = ALPHA1, ALPHA2 = ALPHA2, OMEGL = OMEGL,
+                             GG = GG, BETA = BETA, TAU5= TAU5, 
+                             RHOG=RHOG, material=material, 
+                             min_wvl=min_wvl, max_wvl=max_wvl)
+    
+
+    dni_spectrum = spectral_property(smarts_res[smarts_res.keys()[1]],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    dhi_spectrum = spectral_property(smarts_res[smarts_res.keys()[2]],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    ghi_spectrum = spectral_property(smarts_res[smarts_res.keys()[3]],
+                                     smarts_res['Wvlgth'], interpolation='linear')
+    
+    return (dni_spectrum, dhi_spectrum, ghi_spectrum)
+
+
+
+def spectral_albedo_smarts_SRRL(YEAR, MONTH, DAY, HOUR, ZONE,
+                                LATIT, LONGIT, ALTIT,
+                                RH, TAIR, SEASON, TDAY, SPR, W,
+                                TILT, WAZIM, HEIGHT, 
+                                ALPHA1, ALPHA2, OMEGL, GG, BETA, TAU5,
+                                RHOG, material,
+                                IOUT='30 31', min_wvl='280', max_wvl='4000'):
+ 
+    import pySMARTS
+
+    smarts_res = pySMARTS.SMARTSSRRL(IOUT=IOUT, YEAR=YEAR,MONTH=MONTH,DAY=DAY,HOUR=HOUR, ZONE=ZONE,
+                            LATIT=LATIT, LONGIT=LONGIT, ALTIT=ALTIT, 
+                             RH=RH, TAIR=TAIR, SEASON=SEASON, TDAY=TDAY, SPR=SPR, W=W, 
+                             TILT=TILT, WAZIM=WAZIM, HEIGHT=HEIGHT,
+                             ALPHA1 = ALPHA1, ALPHA2 = ALPHA2, OMEGL = OMEGL,
+                             GG = GG, BETA = BETA, TAU5= TAU5,
+                             RHOG=RHOG, material=material, 
+                             min_wvl=min_wvl, max_wvl=max_wvl)
+    
+    return spectral_property(smarts_res['Zonal_ground_reflectance'],
+                             smarts_res['Wvlgth'], interpolation='linear')
+    
+
+def gen_spectra(idx=None, metdata=None, spectra_folder=None, scale_spectra=False,
+                scale_albedo=False, scale_albedo_nonspectral_sim=False ):
+    
+    # Extract data from metdata
+    dni = metdata.dni[idx]
+    dhi = metdata.dhi[idx]
+    ghi = metdata.ghi[idx]
+    alb = metdata.albedo[idx]
+    solpos = metdata.solpos.iloc[idx]
+    zen = float(solpos.zenith)
+    azm = float(solpos.azimuth) - 180
+    #lat = metdata.latitude
+    #lon = metdata.longitude
+    #elev = metdata.elevation / 1000
+    #t = metdata.datetime[idx]
+
+    # Verify sun up
+    if zen > 90:
+        print("Sun below horizon. Skipping.")
+        return None
+    
+    # Define file suffix
+    # -- CHANGE --
+    suffix = f'_{idx:04}.txt'
+    
+    # Generate/Load dni and dhi
+    dni_file = os.path.join(spectra_folder, "dni"+suffix)
+    dhi_file = os.path.join(spectra_folder, "dhi"+suffix)
+    ghi_file = os.path.join(spectra_folder, "ghi"+suffix)
+    spectral_dni, spectral_dhi, spectral_ghi = spectral_irradiance_smarts(zen, azm, min_wavelength=280)
+    
+    # If specifed, scale the irradiance spectra based on their respective
+    # measured value.
+    if scale_spectra:
+        
+        dni_scale = dni / spectral_dni.df.apply(lambda g: integrate.trapz(spectral_dni.df.value, x=spectral_dni.df.index))
+        dhi_scale = dhi / spectral_dhi.df.apply(lambda g: integrate.trapz(spectral_dhi.df.value, x=spectral_dhi.df.index))
+        ghi_scale = ghi / spectral_ghi.df.apply(lambda g: integrate.trapz(spectral_ghi.df.value, x=spectral_ghi.df.index))
+        # dni_scale = dni / (10*np.sum(spectral_dni[range(280, 4000, 10)]))
+        # dhi_scale = dhi / (10*np.sum(spectral_dhi[range(280, 4000, 10)]))
+        # ghi_scale = ghi / (10*np.sum(spectral_ghi[range(280, 2501, 10)]))
+        spectral_dni.scale_values(dni_scale.value)
+        spectral_dhi.scale_values(dhi_scale.value)
+        spectral_ghi.scale_values(ghi_scale.value)
+    
+    # Write irradiance spectra
+    #'''
+    spectral_dni.to_file(dni_file)
+    spectral_dhi.to_file(dhi_file)
+    spectral_ghi.to_file(ghi_file)
+    #'''
+    
+    # Generate/Load albedo
+    alb_file = os.path.join(spectra_folder, "alb"+suffix)
+    spectral_alb = spectral_albedo_smarts(zen, azm, 'Gravel', min_wavelength=280)
+ 
+    # If specifed, scale the spectral albedo to have a mean value matching the
+    # measured albedo.
+    if scale_albedo:
+        # option A
+        denom = spectral_alb.df.value * spectral_ghi.df.value
+        # option B
+        #denom = spectral_alb.df
+
+        # TODO:
+        # Add test to
+        if alb > 1 or alb == 0:
+            print("albedo measured is an incorrect number, not scaling albedo generated")
+        else:
+            alb_scale = alb / denom.apply(lambda g: integrate.trapz(denom.values, x=spectral_alb.df.index))
+            spectral_alb.scale_values(alb_scale.values)
+
+    if scale_albedo_nonspectral_sim:
+        sim_alb = np.sum(spectral_alb[range(280, 2501, 10)] * spectral_ghi[range(280, 2501, 10)])/np.sum(spectral_ghi[range(280, 2501, 10)])
+        
+        if alb > 1:
+            print("albedo measured is an incorrect number, not scaling albedo generated")
+        else:
+            alb_scale = alb / sim_alb
+            spectral_alb.scale_values(alb_scale)
+            print(alb, sim_alb, alb_scale)
+    
+    # Write albedo file
+    spectral_alb.to_file(alb_file)
+    
+    return (spectral_alb, spectral_dni, spectral_dhi)
\ No newline at end of file
diff --git a/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.ipynb b/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.ipynb
index 0d3bbc35..80b114fc 100644
--- a/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.ipynb
+++ b/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.ipynb
@@ -43,7 +43,7 @@
    "source": [
     "### 1. Run an hourly simulation\n",
     "\n",
-    "This will generate the results over which we will perform the mismatch analysis. Here we are doing only 1 day to make this 'fater'."
+    "This will generate the results over which we will perform the mismatch analysis. Here we are doing only 1 day to make this faster."
    ]
   },
   {
diff --git a/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.py b/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.py
index a192006c..d6f69984 100644
--- a/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.py
+++ b/docs/tutorials/8 - Advanced topics - Calculating Power Output and Electrical Mismatch.py
@@ -31,7 +31,7 @@
 
 # ### 1. Run an hourly simulation
 # 
-# This will generate the results over which we will perform the mismatch analysis. Here we are doing only 1 day to make this 'fater'.
+# This will generate the results over which we will perform the mismatch analysis. Here we are doing only 1 day to make this faster.
 
 # In[6]:
 
diff --git a/requirements.txt b/requirements.txt
index 177a4b58..9b1001ac 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -17,6 +17,7 @@ pvlib==0.7.1
 pvmismatch==4.1
 py==1.10.0
 pyparsing==2.4.7
+pysmarts==0.0.1
 pytest==5.4.1
 pytest-cov==2.8.1
 python-dateutil==2.8.1
diff --git a/setup.py b/setup.py
index 998143d7..a20b82cb 100644
--- a/setup.py
+++ b/setup.py
@@ -120,7 +120,12 @@
     # installed, specify them here.  If using Python 2.6 or less, then these
     # have to be included in MANIFEST.in as well.
     package_data={
-        'bifacial_radiance': ['data/ground.rad','data/gencumulativesky.exe','data/module.json','data/default.ini'],
+        'bifacial_radiance': [
+            'data/ground.rad',
+            'data/gencumulativesky.exe',
+            'data/module.json',
+            'data/default.ini',
+            'images/*'],
     },
     include_package_data=True,
     # Although 'package_data' is the preferred approach, in some case you may


Test patch:
diff --git a/tests/test_bifacial_radiance.py b/tests/test_bifacial_radiance.py
index 91352221..9ca0eb2e 100644
--- a/tests/test_bifacial_radiance.py
+++ b/tests/test_bifacial_radiance.py
@@ -198,8 +198,18 @@ def test_1axis_gencumSky():
     demo.exportTrackerDict(trackerdict, savefile = 'results\exportedTrackerDict')
     assert trackerdict[-5.0]['radfile'][0:7] == 'objects' 
     #assert trackerdict[-5.0]['radfile'] == 'objects\\1axis-5.0_1.825_11.42_5.0_10x3_origin0,0.rad'
-#    trackerdict = demo.makeOct1axis(trackerdict=trackerdict) # just run this for one timestep: Jan 1 11am
-#    trackerdict = demo.analysis1axis(trackerdict=trackerdict, modWanted=7, rowWanted=3, sensorsy=2) 
+    minitrackerdict = {}
+    minitrackerdict[list(trackerdict)[0]] = trackerdict[list(trackerdict.keys())[0]]
+    trackerdict = demo.makeOct1axis(trackerdict=minitrackerdict) # just run this for one timestep: Jan 1 11am
+    trackerdict = demo.analysis1axis(trackerdict=trackerdict, modWanted=7, rowWanted=3, sensorsy=2) 
+    assert trackerdict[-5.0]['AnalysisObj'].x[0] == -10.76304
+    modscanfront = {}
+    modscanfront = {'xstart': -5}
+    trackerdict = demo.analysis1axis(trackerdict=trackerdict, modWanted=7, rowWanted=3, sensorsy=2, modscanfront=modscanfront ) 
+    assert trackerdict[-5.0]['AnalysisObj'].x[0] == -5
+
+
+
 
 def test_SceneObj_makeSceneNxR_lowtilt():
     # test _makeSceneNxR(tilt, height, pitch, azimuth = 180, nMods = 20, nRows = 7, radname = None)
@@ -364,4 +374,25 @@ def test_left_label_metdata():
     demo = bifacial_radiance.RadianceObj('test')
     metdata2 = demo.readEPW(epwfile=MET_FILENAME, label='right' )
     pd.testing.assert_frame_equal(metdata1.solpos, metdata2.solpos)
-    assert metdata2.solpos.index[7] == pd.to_datetime('2001-01-01 07:42:00 -7')
\ No newline at end of file
+    assert metdata2.solpos.index[7] == pd.to_datetime('2001-01-01 07:42:00 -7')
+    
+def test_addMaterialGroundRad():  
+    # test set1axis.  requires metdata for boulder. 
+    name = "_test_addMaterial"
+    demo = bifacial_radiance.RadianceObj(name)
+    demo.setGround(0.2)
+    material = 'demoMat'
+    com = "a demonstration material"
+    Rrefl = 0.9
+    Grefl = 0.2
+    Brefl = 0.9
+    demo.addMaterial(material=material, Rrefl=Rrefl, Grefl=Grefl, Brefl=Brefl, comment=com)
+    demo.setGround('demoMat')
+    assert list(demo.ground.Rrefl) == [0.9]
+    Rrefl = 0.45
+    demo.addMaterial(material=material, Rrefl=Rrefl, Grefl=Grefl, Brefl=Brefl, comment=com, rewrite=False)
+    demo.setGround('demoMat')
+    assert list(demo.ground.Rrefl) == [0.9]
+    demo.addMaterial(material=material, Rrefl=Rrefl, Grefl=Grefl, Brefl=Brefl, comment=com, rewrite=True)
+    demo.setGround('demoMat')
+    assert list(demo.ground.Rrefl) == [0.45]
\ No newline at end of file