TJ clear-sky model implementations use sin instead of sind

[kandersolar]

I think all three implementations of the TJ clear-sky model in this repo are flawed in that the arguments to sin are interpreted as radians but should be interpreted as degrees. Here are the relevant lines of code:

Engerer2-separation-model/Engerer2Separation.py

Lines 128 to 130 in 8749240

	a_tj = 1160 + 75 * np.sin((360 * (day_of_year - 275)) / 365)
	k_tj = 0.174 + 0.035 * np.sin((360 * (day_of_year - 100)) / 365)
	c_tj = 0.095 + 0.04 * np.sin((360 * (day_of_year - 100)) / 365)

Engerer2-separation-model/Engerer2Separation.m

Lines 122 to 124 in 8749240

	A = 1160 + 75 .* sin((360.*(doy-275))./365);
	k = 0.174 + 0.035 .* sin((360.*(doy-100))./365);
	C = 0.095 + 0.04 .* sin((360.*(doy-100))./365);

Engerer2-separation-model/Engerer2Separation.R

Lines 97 to 99 in 8749240

	A = 1160 + 75 * sin((360 * (doy - 275)) / 365)
	k = 0.174 + 0.035 * sin((360 * (doy - 100)) / 365)
	C = 0.095 + 0.04 * sin((360 * (doy - 100)) / 365)

As an aside, I think the implementation in the clear-sky-models repo has the same problem: https://github.com/JamieMBright/clear-sky-models/blob/86fe607e9ba66226be6102895491b7d756473a7e/R/1-TJ.R#L26-L28

The 360/365 coefficient seems like it is intended to produce one period per year in degrees. Treating it as radians introduces an undesirable amplitude variation of period ~6 days, as these two screenshots show:

Changing the 360s to 2*np.pi seems to resolve the issue in the python implementation. I did not check the R or MATLAB versions.

As far as I can tell, these equations for A, k, and C do not come from the 1957 Threlkeld & Jordan paper but rather from G Masters' "Renewable and Efficient Electric Power Systems", which contains this example (Example 7.8). It only gives the correct answer if the argument to sin is interpreted as degrees:

[kandersolar]

Perhaps worth pointing out that the parameterization script has the same:

Engerer2-separation-model/ParameteriseEngerer2Separation.m

Lines 141 to 143 in d09b69d

	A = 1160 + 75 .* sin((360.*(doy-275))./365);
	k = 0.174 + 0.035 .* sin((360.*(doy-100))./365);
	C = 0.095 + 0.04 .* sin((360.*(doy-100))./365);

[JamieMBright]

Yeah you are of course right.
Thanks for finding this. I shall push a change now with the relevant updates.

[JamieMBright]

I have committed the changes to master.
Thanks for the open collaborative spirit here.
It might well be updating this to have a more consistent clear-sky model worldwide following a model from here: https://doi.org/10.1016/j.rser.2020.110087

[kandersolar]

Thanks @JamieMBright!

I think there is unfortunately no such function np.sind (or at least my numpy doesn't have it). So perhaps np.sin(np.radians(...)) or np.sin(np.pi/180 * ...)) should be used in the python example.

What would you suggest about using the parameter values published in https://aip.scitation.org/doi/10.1063/1.5097014? Would you be worried about error because of this issue in how the parameter values were optimized?

[JamieMBright]

Right you are.

Yeah to be honest it would be uncertain now. Changing this code and then using the same parameters could/likely will lead to adverse affects. It won't be catastrophic as you demonstrate the values produced simply have a subtle oscillation rather than anything wild. I would need to revisit/update this work, though the process/method is still identical. Sadly I don't have the data in as clean a format as I once had at my previous job.

I believe sind is fine for both R and matlab. Changing the Python implementation to np.sin(np.deg2rad(...))

[AdamRJensen]

@JamieMBright Do you know if the coefficients in the 2015 paper were based on the same incorrect implementation of the clear-sky model?

May I also suggest that a message be placed in the readme file of the repository letting users know of this known issue / change?