If there is only one observation in the day or if \(\small t\) falls between sunrise and the first observation or between the last observation and sunset then we just use the ratio for the nearest observation (\(\small R_1\)):\(\small InstPAR(t) = R_1*\sin(elevation(t))\). Below I dissect the functional form of \(\small \sin(elevation(t))\) based off NOAA's calculation of solar elevation\(^2\). $$\small \sin(elevation(t))=\sin(90-\arccos(\sin(Lat)*\sin(SunDecline)+\cos(Lat)*\cos(Sundecline)*\cos(t*360+EqOfTime/4+Lon+180)))$$ $$\small \sin(elevation(t))=\cos(Arccos(\sin(Lat)*\sin(SunDecline)+\cos(Lat)*\cos(Sundecline)*-\cos(t*360+EqOfTime/4+Lon)))$$ $$\small \sin(elevation(t))=\sin(Lat)*\sin(SunDecline)-\cos(Lat)*\cos(Sundecline)*\cos(t*360+EqOfTime/4+Lon)$$ Since SunDecline and EqOfTime are approximately constant during a given day and since Lat and Lon don't change at a location, if we set \(\small \sin(Lat)*\sin(SunDecline)=A, \hspace{4pt} \cos(Lat)*\cos(Sundecline)=B,\hspace{4pt}\) and \(\small EqOfTime/4+Lon=C\) we get \(\hspace{4pt} \small \sin(elevation(t))=A-B*\cos(t*360+C)\) where \(\small A,B,C\) are day and location specific constants. So, substituting that into our \(\small InstPAR(t)\) function we get: $$\small InstPAR(t) = (R_2 * \frac{t-T_1}{T_2-T_1}+R_1*\frac{T_2-t}{T_2-T_1})*(A-B*\cos(t*360+C))$$ when there are daylight observations on either side of \(\small t\). When there is only one daylight observation close to \(\small t\) then: $$\small InstPAR(t) = R_1*(A-B*\cos(t*360+C))$$

So, to calculate PAR for a given day at a specified location with \(\small n\) daylight observations we integrate \(\small InstPAR(t)\) from sunrise to sunset: $$\small \int_{Sunrise}^{T_1} R_1*(A-B*\cos(t*360+C))dt + \int_{T_1}^{T_2} (R_2 * \frac{t-T_1}{T_2-T_1}+R_1*\frac{T_2-t}{T_2-T_1})*(A-B*\cos(t*360+C))dt +...$$ $$\small + \int_{T_{n-1}}^{T_n} (R_n *\frac{t-T_{n-1}}{T_n-T_{n-1}}+R_{n-1}*\frac{T_n-t}{T_n-T_{n-1}})*(A-B*\cos(t*360+C))dt + \int_{T_n}^{Sunset} R_n*(A-B*\cos(t*360+C))dt$$

2. https://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html

3. Liang S., Zhang X., Xiao Z., Cheng J., Liu Q., Zhao X. (2014) Incident Photosynthetic Active Radiation. In: Global LAnd Surface Satellite (GLASS) Products. SpringerBriefs in Earth Sciences. Springer, Cham

Input values except for time are self explanatory. The time variable represents the time in GMT standardized such that a full, 24 hour, day is equal to length 1. In the data used, Incident PAR values were collected every 3 hours for a total of 8 times per day. As such, you can set time equal to any time incident PAR would be collected. The time variable is only used in the first plot to depict the elevation needed by the Wang et al. method for the ratio of interpolated values to exhibit a linear relationship between observations. In plots 2 and 3, for each observation at a time where there would be sunlight, a number between 0 and 1 is randomly selected to be the PAR ratio for that observation. That number is then multiplied by the theoretical maximum PAR value at that solar elevation to create the PAR value at that observation. Note: random gerenation is used here as these plots are just for visualization puropses for how the interpolation methods work and differ. You will notice that Wang et al.'s algorithm works great for the vast majority of observations. However, the closer we get to the poles the more inaccurate Wang's algorithms get (try moving latitude to -77 or so with the rest of the options in their defult settings).