Merge pull request #105 from cse-sim/solarinterp

Solarinterp

Author: chipbarnaby

src/CGCOMP.CPP

@@ -508,13 +508,10 @@ const double G0STDX = g0Std;	// CSE accel. of gravity, ft/sec2
 // class ANDAT -- airnet data
 void ANDAT::ad_ClearResults()
 {
-	ad_delP = 0.;
+	ad_delP = ad_dmdp = 0.;
 	ad_mdotP = ad_mdotB = ad_mdotX = 0.;
-	ad_dmdp = 0.;
-	ad_xDelpF = 0.;
-	ad_xMbm = 0.;
-	ad_tdFan = 0.;
-	ad_pFan = 0.;
+	// ad_xDelpF = ad_xMbm = 0.;	do not clear (values retained throughout hour)
+	ad_tdFan = ad_pFan = 0.;
 
 }	// ANDAT::ad_ClearResults
 //-----------------------------------------------------------------------------

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

@@ -81,6 +81,14 @@ RC TOPRAT::tp_WthrInit()  	// Initialize weather data/ display any error message
 	return rc;
 }		// TOPRAT::tp_WthrInit
 //---------------------------------------------------------------------------
+RC TOPRAT::tp_WthrBegDay()
+{
+	RC rc = RCOK;
+	rc = Wfile.wf_UpSampleDay(jDay);
+
+	return rc;
+}		// TOPRAT::tp_WthrBegDay
+//---------------------------------------------------------------------------
 RC TOPRAT::tp_WthrBegHour()		// start-hour weather stuff: read file, set up public Wthr and Top members.
 // tp_WthrBegSubhr must also be called, after this function.
 
@@ -259,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
@@ -396,31 +414,46 @@ 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)
+#if defined( _DEBUG)
+	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");
+#endif
+#else
+	if (!rc)			// if wfRead failed, eg date not in file or type without design days
+		wd_Adjust(iHrST);
+#endif
+
+	return rc;
+}		// WDHR::wd_WfReader
+//-----------------------------------------------------------------------------
+void WDHR::wd_Adjust(		// apply adjustments to weather data
+	int iHrST)		// standard time hour of day (0-23)
+{
 
-// 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)
+	// 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. */
 
-	// 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);
+	// 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?
 
+	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
+}		// 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/CNGUTS.CPP

@@ -504,6 +504,7 @@ LOCAL RC FC NEAR doBegIvl()	// simulation run start-of-interval processing: init
 		Top.tp_DoDateDowStuff();
 		// CAUTION: Top.iHr may be --'d or ++'d in tp_DoDTStuff in C_IVLCH_H case if DT starts or ends.
 		Top.isSolarCalcDay = Top.isBegMonth;     		// do solar calcs on warmup/run/mon 1st day and dsn day 1st rep.
+		Top.tp_WthrBegDay();
 		cgReportsDaySetup();		// init zones & Top re date-dependent reports and exports.
 				     				// set zp-> lists; may print heads.  Uses jDay. cgresult.cpp.
 

src/CNRECS.DEF

@@ -517,20 +517,17 @@ RECORD TOPRAT "top" *RAT	/* top level RAT: contains control info and all once-on
     *s POWER radBeamShAv	// .. current beam subhour average power, interpolated, Btuh/ft2
   #ifdef SOLAVNEND // cndefns.h (via cnglob.h) 1-18-94: only if computing & using end-ivl as well as ivl avg solar values
   o  *s POWER radBeamSh		// .. end subhour instantaneous value, interpolated
-  #else
-    FLOAT radBeamShSpare
   #endif
     *h POWER radDiffHrAv 	// diffuse irradiance on horizontal surface, hour energy = average power, from weather file
     *h POWER radDiffPvHrAv	// .. flavors as for radBeam
     *h POWER radDiffNxHrAv	// ..
     *s POWER radDiffShAv	// .. current diffuse subhour power, interpolated by cgwthr.cpp, Btuh/ft2
   #ifdef SOLAVNEND // cndefns.h
   o  *s POWER radDiffSh		// .. end subhour instantaneous value, interpolated
-  #else
-    FLOAT radDiffShSpare
   #endif
 
     // other weather: weather file contains instantaneous values for END OF HOUR; program needs only end-interval values.
+	*declare "RC tp_WthrBegDay();"
     *declare "RC tp_WthrBegHour();"
     *declare "RC tp_WthrBegSubhr();"
     *h TEMP tDbOHr			// outdoor dry bulb temp at end of hour, from wthr file, deg F.
@@ -690,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);"
@@ -769,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 */
   {