Solar calc reorg WIP; test ref files updated

Author: chipbarnaby

src/CGSOLAR.CPP

@@ -418,8 +418,9 @@ void FC makHrSgt(				// make solar tables for an hour for current month
 	// Re-initialize slpak for current day as assumed by cgwfread() called from cnguts.cpp each hour b4 this fcn.
 	// Don't know a) if other parts of program depend on it (else do in wswfread each call)
 	//            b) if would be faster to push & pop the data it sets.
-	slday( Top.jDayST,		// julian date of simulation (standard time), cnguts.cpp.  Use of ST matches call in cgwthr.cpp.
-		   SLTMLST );		// say use local standard time
+	slday(Top.jDayST,		// julian date of simulation (standard time), cnguts.cpp.  Use of ST matches call in cgwthr.cpp.
+			SLTMLST,		// say use local standard time
+			1);				// skip if no day change
 		// slpak.cpp: set up curr SLLOCDAT struct (which is cse.cpp:Locsolar) for given day:
 		// sets info re declination of earth's axis, hourly sunupf[], dircos[], slazm[], etc.
 		// Used by most other slpak calls.
@@ -873,7 +874,7 @@ void SgThruWin::tw_Doit()
 
 		// target gain for SGDISTs
 
-		for (sgi = 0; sgi < tw_xr->x.nsgdist; sgi++)	// explicit user-entered targets
+		for (int sgi = 0; sgi < tw_xr->x.nsgdist; sgi++)	// explicit user-entered targets
 		{
 			SGDIST* sgd = tw_xr->x.sgdist + sgi; 		// point to sgdist
 			if (sgd->sd_targTy == SGDTTSURFI || sgd->sd_targTy == SGDTTSURFO)

src/CGWTHR.CPP

@@ -84,7 +84,7 @@ RC TOPRAT::tp_WthrInit()  	// Initialize weather data/ display any error message
 RC TOPRAT::tp_WthrBegDay()
 {
 	RC rc = RCOK;
-
+	rc = Wfile.wf_UpSampleDay(jDay);
 
 	return rc;
 }		// TOPRAT::tp_WthrBegDay
@@ -267,6 +267,16 @@ o   else { aes = 0.f; bes = 1.5f - g; ces = g - .5f; }
 	radBeamShAv = (a * radBeamPvHrAv + b * radBeamHrAv + c * radBeamNxHrAv) * beamAdj;
 	radDiffShAv = (a * radDiffPvHrAv + b * radDiffHrAv + c * radDiffNxHrAv) * diffAdj;
 
+#if 0
+	float radBeamX, radDiffX;
+	Wfile.wf_GetSubhrSolar(iHr, iSubhr, radBeamX, radDiffX);
+	if (frDiff(radBeamShAv, radBeamX) > .001f || frDiff(radDiffShAv, radDiffX) > .001f)
+	{
+		printf("\nSolar mismatch");
+		Wfile.wf_UpSampleDay(jDay);
+	}
+#endif
+
 // Interpolate end-subhour instantaneous values for same data
 
 #ifdef SOLAVNEND // undef in cndefns.h
@@ -404,31 +414,44 @@ RC WDHR::wd_WfReader( 	// read an hour's weather data and make adjustments
 		rc = pWF->wf_Read( this, jDayST, iHrST, WRN|wfOp);	// Read hour's data from weather file
 		break;
 	}
-	if (rc)			// if wfRead failed, eg date not in file or type without design days
-		return rc;
-
-// Optional anisotropic sky adjustment
-// Do NOT adjust wd_DNI, wd_DHI
-	if (Top.skyModel==C_SKYMODCH_ANISO)
-		slaniso( &wd_bmrad, &wd_dfrad, iHrST);	// adjust data in place to approx model sky brighter
-			       								//   near sun, etc, while still computing with sun's
-												//	 beam and isotropic diffuse radiation. slpak.cpp. */
+#if defined( SOLARFIX)
+	float bmrad = wd_DNI;
+	float dfrad = wd_DHI;
+	if (Top.skyModel == C_SKYMODCH_ANISO)
+		slaniso(&bmrad, &dfrad,iHrST);
+	if (frDiff(bmrad, wd_bmrad) > .0001f
+	 || frDiff(dfrad, wd_dfrad) > 0.0001f)
+		printf("\nDiff");
+#else
+	if (!rc)			// if wfRead failed, eg date not in file or type without design days
+		wd_Adjust(iHrST);
+#endif
 
-// other adjustments
-	wd_bmrad *= Top.radBeamF; 	// Adjust radiation values for possible user input factors, defaults 1.0.
-	wd_dfrad *= Top.radDiffF; 	// ..
-	wd_wndSpd = Top.windF * max( wd_wndSpd, Top.windSpeedMin);
-					// apply wind speed min (user input, default = .5)
-					//   and adjust by factor (per user input, default .25)
+	return rc;
+}		// WDHR::wd_WfReader
+//-----------------------------------------------------------------------------
+void WDHR::wd_Adjust(		// apply adjustments to weather data
+	int iHrST)		// standard time hour of day (0-23)
+{
 
-	// Ensure solar data values are not negative (as found in some CSW files)
-	wd_bmrad = max(wd_bmrad, 0.f);
-	wd_dfrad = max(wd_dfrad, 0.f);
-	wd_glrad = max(wd_glrad, 0.f);
+	// Optional anisotropic sky adjustment
+	// Do NOT adjust wd_DNI, wd_DHI
+	if (Top.skyModel == C_SKYMODCH_ANISO)
+		slaniso(&wd_bmrad, &wd_dfrad, iHrST);	// adjust data in place to approx model sky brighter
+												//   near sun, etc, while still computing with sun's
+												//	 beam and isotropic diffuse radiation. slpak.cpp. */
 
+	// other solar adjustments
+	//   ensure solar data values are not negative (as found in some CSW files)
+	wd_bmrad = max( Top.radBeamF*wd_bmrad, 0.f); 	// Adjust radiation values by possible user input factors, defaults 1.0.
+	wd_dfrad = max( Top.radDiffF*wd_dfrad, 0.f); 	// ..
+	wd_glrad = max(wd_glrad, 0.f);		// TODO: re-derive from adjusted wd_bmrad and wd_dfrad?
 
-	return rc;
-}		// WDHR::wd_WfReader
+	wd_wndSpd = Top.windF * max(wd_wndSpd, Top.windSpeedMin);
+	// apply wind speed min (user input, default = .5)
+	//   and adjust by factor (per user input, default .25)
+
+}		// WDHR::wd_Adjust
 //---------------------------------------------------------------------------
 void FC cgWfDone()   		// Close, clean up hourly simulator weather file
 {

src/CNDEFNS.H

@@ -132,9 +132,9 @@
 #define DIM_POLYGONXYZ 37		// input arrays dimension
 								//   = MAX_POLYLPVERTICIES*3 + 1
 
-#undef SOLARFIX				// define to correct timing of slday() calls
+#define SOLARFIX			// define to correct timing of slday() calls
 							//    and definition of HA, 10-25-2017
-
+							//    plus add'l changes, 8-2020
 #endif	// ifndef _CNEDEFNS_H
 
 // cndefns.h end

src/CNRECS.DEF

@@ -687,6 +687,8 @@ RECORD WFILE "weatherFile" *RAT	// weather file info, one static instance "Wfile
 	*declare "RC wf_EtDecodeHdr( char* hdr, int erOp, float* clrnss, float* turbid, float* atmois);"
 	*declare "RC wf_DecodeHdrFields( char* hdr, struct WFHTAB* wfht0, int erOp);"
 	*declare "RC wf_FillWDYEAR( int erOp=WRN);"
+	*declare "RC wf_UpSampleDay( int jDay);"
+	*declare "void wf_GetSubhrSolar( int iHr, int iSh, float& radBeam, float& radDiff) const;"
 	*declare "RC wf_Read( WDHR* pWd, int jDay, int iHr, int erOp = WRN);"
 	*declare "USI* wf_PackedHrRead( int jDay, int iHr, int erOp = WRN);"
 	*declare "LI wf_PackedHrOffset( int jDay, int iHr, int erOp = WRN);"
@@ -766,6 +768,7 @@ RECORD WDHR "wfdata sub" *SUBSTRUCT	// hourly data substructure for WFDATA
     *declare "void wd_Init( int options=0);"
     *declare "WDHR& Copy( const WDHR& wd, int options=0);"
 	*declare "RC wd_WfReader( BOO nextHour, WFILE* pWF);"
+	*declare "void wd_Adjust( int iHrST);"
 	*declare "RC wd_Unpack( int iH, USI* pHour, int wFileFormat=ET1);"
     *declare "RC wd_EstimateMissingET1( int iHr);"
     *declare "float wd_CalcSkyTemp( int skyModelLW, int iHr);"

src/SLPAK.CPP

@@ -352,10 +352,6 @@ void FC slday( // set up daily data in current SLLOCDAT
   if ((options & 1) && slloccur->doy == doy && slloccur->timetype == timetype)
     return; // already set for doy
 
-#if 0 && defined(_DEBUG)
-	printf("\nslday %d %d", doy, timetype);
-#endif
-
   slloccur->doy = doy;
 
   float dec;
@@ -746,35 +742,35 @@ void FC slaltazm(/* Calculate solar altitude and azimuth for hour */
 } /* slaltazm */
 //======================================================================
 void FC slaniso(/* Adjust beam and diffuse radiation values
-       according to Hay anisotropic sky model */
-
-                float *pbeam, /* Pointer to beam value (Btu/sf): used/replaced. */
-                float *pdiff, /* Pointer to diffuse value (Btu/sf): used/replaced.*/
-                SI ihr)       /* Hour of day (0 - 23, 0 = midnight to 1 AM).
-                           Meaning of hour (LST, solar time, etc.)
-                           depends on how SLLOCDAT was set up by slday(). */
-
-/* must call slday() for day first (uses slloccur members) */
-
-/* Adjusted values of beam and diffuse *REPLACE* the current values */
-
-/* Story: Code assumes sky is modelled as istropic (uniform) hemisphere of
-   sky radiation plus direct sun beam.  But actually, sky is brighter
-   near sun, near horizon, etc.  Hay model approximates reality better
-   with old code by increasing beam to approximate some of the extra sky
-   brightness near sun.  Rob per Chip, 12-89. */
-
-/* Recoded 2-10-89 based on cp4sim equivalent function.  Features --
-  1.  Doesn't bother unless there is a little beam
-  2.  Constrains fb to be .8 max.  This prevents wild values from
-      early morning observations.
-  3.  Derives fd (diffuse factor) from fb rather than from basic
-      Hay formula, giving same result unless fb hit .8 limit.  In all
-      cases, total horiz is preserved */
+	   according to Hay anisotropic sky model */
+
+	float *pbeam, /* Pointer to beam value (Btu/sf): used/replaced. */
+	float *pdiff, /* Pointer to diffuse value (Btu/sf): used/replaced.*/
+	SI ihr)       /* Hour of day (0 - 23, 0 = midnight to 1 AM).
+			   Meaning of hour (LST, solar time, etc.)
+			   depends on how SLLOCDAT was set up by slday(). */
+
+			   /* must call slday() for day first (uses slloccur members) */
+
+			   /* Adjusted values of beam and diffuse *REPLACE* the current values */
+
+			   /* Story: Code assumes sky is modelled as istropic (uniform) hemisphere of
+				  sky radiation plus direct sun beam.  But actually, sky is brighter
+				  near sun, near horizon, etc.  Hay model approximates reality better
+				  with old code by increasing beam to approximate some of the extra sky
+				  brightness near sun.  Rob per Chip, 12-89. */
+
+				  /* Recoded 2-10-89 based on cp4sim equivalent function.  Features --
+					1.  Doesn't bother unless there is a little beam
+					2.  Constrains fb to be .8 max.  This prevents wild values from
+						early morning observations.
+					3.  Derives fd (diffuse factor) from fb rather than from basic
+						Hay formula, giving same result unless fb hit .8 limit.  In all
+						cases, total horiz is preserved */
 {
-  SI pos1, pos2;
-  SI ihx;
-  float c1, c2, cosi, f, fb, fd;
+	SI pos1, pos2;
+	SI ihx;
+	float c1, c2, cosi, f, fb, fd;
 
   if (*pbeam > 5.f) /* if a little beam */
   {

src/WFPAK.CPP

@@ -12,6 +12,7 @@
 #include "ancrec.h"	// record: base class for rccn.h classes
 #include "rccn.h"	// TOPRAT WFILE WFDATA
 #include "psychro.h"
+#include "slpak.h"
 
 #include "yacam.h"	// class YACAM
 #include "tdpak.h"
@@ -445,8 +446,10 @@ struct WDSUBHRFACTORS		// constant factors for weather data interpolation
 class WDSUBHR		// derived subhourly weather data
 {
 public:
-	float wsh_radBeam;
-	float wsh_radDiff;
+	void wsh_ClearSolar()
+	{	wsh_radBeamAv = wsh_radDiffAv = 0;	}
+	float wsh_radBeamAv;
+	float wsh_radDiffAv;
 };		// class WDSUBHR
 //----------------------------------------------------------------------------
 struct VHR		// sort helper for wdy_Stats
@@ -501,18 +504,19 @@ public:
 	float wdy_TaDbAvg( int jDay1, int jDay2);
 	RC wdy_UpSampleSetup(int nSh);
 	RC wdy_UpSampleDay(int jDay, int nSh);
-	RC wdy_UpSampleHr(int jDay, int nSh);
+	RC wdy_UpSampleHr(WDSUBHR wdsh[], int jDay, int iHr, int nSh);
 };		// class WDYEAR
 //-----------------------------------------------------------------------------
 WDYEAR::WDYEAR()			// c'tor
-	: wdy_wShF(NULL)
+	: wdy_wShF(NULL), wdy_shDay(NULL)
 {
 	wdy_Init();
 }		// WDYEAR::WDYEAR
 //-----------------------------------------------------------------------------
 WDYEAR::~WDYEAR()			// d'tor
 {
 	delete[] wdy_wShF;
+	delete[] wdy_shDay;
 }		// WDYEAR::~WDYEAR
 //-----------------------------------------------------------------------------
 void WDYEAR::wdy_Init(
@@ -589,10 +593,18 @@ RC WDYEAR::wdy_Fill(	// read weather data for entire file; compute averages etc.
 	{	int monLen = monLens[ iMon]  + (wdy_isLeap && iMon==2);
 		for (iDay=0; iDay<monLen; iDay++)
 		{	jDay++;
+#if 1
+			slday(jDay, SLTMLST);
+#endif
 			for (int iHr=0; iHr<24; iHr++)
-			{	rc |= pWF->wf_Read( &wdy_Hr( jDay, iHr), jDay, iHr, erOp|WF_FORCEREAD);
+			{	WDHR* pWd = &wdy_Hr(jDay, iHr);
+				rc |= pWF->wf_Read( pWd, jDay, iHr, erOp|WF_FORCEREAD);
 				if (rc)
 					break;
+#if defined( SOLARFIX)
+				pWd->wd_Adjust( iHr);	// apply adjusments per user input
+										//   aniso, solar factors, wind factor
+#endif
 			}
 			if (rc)
 				break;
@@ -790,6 +802,8 @@ RC WDYEAR::wdy_UpSampleSetup(
 
 	delete[] wdy_wShF;
 	wdy_wShF = new WDSUBHRFACTORS[nSh];
+	delete[] wdy_shDay;
+	wdy_shDay = new WDSUBHR[ 24*nSh];
 
 	for (int iSh = 0; iSh < nSh; iSh++)
 		wdy_wShF[iSh].wds_Setup(iSh, nSh);
@@ -798,18 +812,39 @@ RC WDYEAR::wdy_UpSampleSetup(
 
 }		// WDYEAR::wdy_UpSampleSetup
 //-----------------------------------------------------------------------------
+RC WFILE::wf_UpSampleDay(
+	int jDay)
+{
+	return wf_pWDY ? wf_pWDY->wdy_UpSampleDay(jDay, Top.nSubSteps)
+		           : RCBAD;
+}		// WFILE::wf_UpSampleDay
+//-----------------------------------------------------------------------------
+void WFILE::wf_GetSubhrSolar(
+	int iHr,
+	int iSh,
+	float& radBeam,
+	float& radDiff) const
+{
+	const WDSUBHR* pWDSH = wf_pWDY->wdy_shDay + iHr*Top.nSubSteps + iSh;
+	radBeam = pWDSH->wsh_radBeamAv;
+	radDiff = pWDSH->wsh_radDiffAv;
+}
+//-----------------------------------------------------------------------------
 RC WDYEAR::wdy_UpSampleDay(
 	int jDay,	// day of year (1-365/366)
 	int nSh)
 {
-
+	RC rc = RCOK;
 	for (int iHr = 0; iHr < 24; iHr++)
 	{
+		WDSUBHR* pWDSH = wdy_shDay + iHr * nSh;
+		rc |= wdy_UpSampleHr(pWDSH, jDay, iHr, nSh);
 	}
+	return rc;
 }		// WDYEAR::wdy_UpSampleDay
 //-----------------------------------------------------------------------------
 RC WDYEAR::wdy_UpSampleHr(
-	XXX wdSh,
+	WDSUBHR wdSh[],	// returned: subhour values
 	int jDay,		// day of year (1-365/366)
 	int iHr,		// hour of day, 0-23
 	int nSh)		// # of subhrs per hour
@@ -827,8 +862,10 @@ RC WDYEAR::wdy_UpSampleHr(
 	}
 
 	if (radBeam[1] < .1f && radDiff[1] < .1f)
-	{	// substep 0
-		continue;
+	{
+		for (int iSh = 0; iSh < nSh; iSh++)
+			wdSh[iSh].wsh_ClearSolar();
+		return rc;
 	}
 
 	// Compute adjustment factors (constant over hour) to make subhour values add up to given hour value.
@@ -871,8 +908,9 @@ RC WDYEAR::wdy_UpSampleHr(
 
 		// Interpolate subhour average power values for data with average values in file.
 		//   Multiply by subhrDur for energy.
-		wdSh[ iSh].radBeamAv = VIProd(radBeam, 3, wf[iSh].slrF) * beamAdj;
-		wdSh[ iSh].radDiffAv = VIProd(radDiff, 3, wf[iSh].slrF) * diffAdj;
+		const float* slrF = wdy_wShF[iSh].slrF;
+		wdSh[ iSh].wsh_radBeamAv = VIProd( radBeam, 3, slrF) * beamAdj;
+		wdSh[ iSh].wsh_radDiffAv = VIProd( radDiff, 3, slrF) * diffAdj;
 
 	}
 
@@ -1278,7 +1316,7 @@ WFILE::WFILE( basAnc *b, TI i, SI noZ /*=0*/)	// c'tor to help insure init is ca
 }		// WFILE::WFILE
 //----------------------------------------------------------------------------
 WFILE::WFILE()			// default c'tor (insurance)
-	: record( NULL, 0, 1), yac( NULL), yacTDV( NULL)
+	: record( NULL, 0, 1), yac( NULL), yacTDV( NULL), wf_pWDY( NULL)
 {	wf_Init();
 }		// WFILE::WFILE
 //----------------------------------------------------------------------------
@@ -2665,7 +2703,7 @@ x			printf( "mismatch\n");
 	return rc;
 }		// WDHR::wd_EPWReadHr
 #undef _C
-
+//-----------------------------------------------------------------------------
 RC WFILE::wf_FixJday(DOY& jDay, int begDay)
 {
 	if (jd1 > 0				// if weather file dates known