Skip to content

Incidence Angle Modifier zero instead of np.nan (#338) #339

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 16 commits into from
Aug 1, 2017

Conversation

jkfm
Copy link
Contributor

@jkfm jkfm commented Jun 14, 2017

This PR changes the three IAM-models to return 0 instead of np.nan for aoi >= 90 (ashraeiam, physicaliam) and <0 (sapm_aoi_loss).

I also updated the tests to account for the changes. The test for physicaliam now shows, that angles of -90 and 90 return 0, whereas an angle of 0 returns 1, in accordance to the model description in De Soto et al..

Comments are welcome!

@adriesse
Copy link
Member

Could we have the identical code in each of the three functions to do the same thing? That is, to deal with the negative values and the larger than 90 angles? I don't know whether np.maximum, np.where or np.clip is better, but I would lean toward clip.

The SAPM model doesn't like negative angles, the others are fine because of the cosine. That should be handled consistently as well.

@wholmgren wholmgren added this to the 0.5.0 milestone Jun 14, 2017
@wholmgren wholmgren added the api label Jun 14, 2017
Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jkfm. Can you also make sure that the tests assert than the functions preserve nan values? I think they currently do, but this will prevent someone from inadvertently adding e.g. fillna to the functions in the future.

@@ -811,6 +812,10 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
spa
ashraeiam
'''
zeroang = 1e-06

aoi = np.where(aoi == 0, zeroang, aoi)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did you add this line? I understand what it does, but not why we need it now when we didn't previously.

It's not clear to me that we actually need zeroang at all. I think its only purpose is to avoid infs, but maybe those be more cleanly handled with np.where or similar.

Copy link
Contributor Author

@jkfm jkfm Jun 16, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I understood it, if the input angle is 0, tau will be nan, which is why we substitute it with zeroang.
As iam is defined as tau(theta) / tau(0°), we need zeroang anyway to calculate tau(0°). That's why I moved the variable definition to the beginning of the function.

I guess we could define tau0 as a constant in the function and handle aoi == 0 with np.where after the calculations are done, but I feel like the current function is more intuitive to understand.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should keep zeroang and for input theta < zeroang, physicaliam should return tau0. Passing theta=0 through will return nan rather than zero. Line 821 computes tau with the leading factor exp(-KL/cos(theta)) which is why we should return tau0 for small theta.

tau0 is computed awkwardly. I'd replace line 829 - 833 with the limit value which we should have implemented in Matlab:
tau0 = exp(-KL)*(1 - ((1-n)/(1+n))^2)
We used zeroang instead in Matlab which is probably why its in pvlib-python

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, now I see that you've changed the test so that an input of 0 returns 1 instead of np.nan. I should have caught this when I wrote the test.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am surprised that the model would not work for an angle of exactly zero degrees. Could this have something to do with the equation discrepancy mentioned in the comments further above?

Copy link
Contributor Author

@jkfm jkfm Jul 3, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it has to do with the function discrepancy. The unmodified function for theta_r is not defined for angles > ~40°.

The tau(theta)-function in turn simply has a discontinuity for theta = 0. It seems this was ignored in the paper...

@@ -732,7 +732,7 @@ def ashraeiam(aoi, b=0.05):

iam = 1 - b*((1/np.cos(np.radians(aoi)) - 1))

iam = np.where(np.abs(aoi) >= 90, np.nan, iam)
iam = np.where(np.abs(aoi) >= 90, 0, iam)
iam = np.maximum(0, iam)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the np.maximum line still necessary?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is, because for angles close to +/-90° the function would return a negative value for the iam (the intercept for b=0.05 is at 87.21°: https://www.desmos.com/calculator/y8wvzrx3b9).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, the np.maximum line is needed

@wholmgren
Copy link
Member

@adriesse brings up a separate but closely related point regarding the meaning of negative values of aoi in pvlib. I tend to agree that the aoi loss functions should all start with something like aoi = np.abs(aoi), but that's open for discussion. I thought that we had previously discussed this elsewhere, but I can't find the discussion. If we're going to change it, now is the time.

@cwhanse
Copy link
Member

cwhanse commented Jun 14, 2017 via email

@jkfm jkfm changed the title Modelchain zero instead of npnan (#338) Incidence Angle Modifier zero instead of np.nan (#338) Jun 16, 2017
Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jkfm.

If we're not going to treat negative aoi as being equivalent to positive aoi, then my preference is for nans for both AOI < 0 and AOI > 90. Another option is 0 for either or both of those ranges.

In any case, we should update the Variables and Symbols page with the allowable AOI domain. http://pvlib-python.readthedocs.io/en/latest/variables_style_rules.html

@@ -811,6 +812,10 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
spa
ashraeiam
'''
zeroang = 1e-06

aoi = np.where(aoi == 0, zeroang, aoi)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, now I see that you've changed the test so that an input of 0 returns 1 instead of np.nan. I should have caught this when I wrote the test.

@cwhanse
Copy link
Member

cwhanse commented Jun 30, 2017 via email

@jkfm
Copy link
Contributor Author

jkfm commented Jul 3, 2017

I would also prefer IAM to return zero for angles outside -90/90 degrees.

@@ -732,7 +732,7 @@ def ashraeiam(aoi, b=0.05):

iam = 1 - b*((1/np.cos(np.radians(aoi)) - 1))

iam = np.where(np.abs(aoi) >= 90, np.nan, iam)
iam = np.where(np.abs(aoi) >= 90, 0, iam)
iam = np.maximum(0, iam)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, the np.maximum line is needed

@@ -811,6 +812,10 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):
spa
ashraeiam
'''
zeroang = 1e-06

aoi = np.where(aoi == 0, zeroang, aoi)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should keep zeroang and for input theta < zeroang, physicaliam should return tau0. Passing theta=0 through will return nan rather than zero. Line 821 computes tau with the leading factor exp(-KL/cos(theta)) which is why we should return tau0 for small theta.

tau0 is computed awkwardly. I'd replace line 829 - 833 with the limit value which we should have implemented in Matlab:
tau0 = exp(-KL)*(1 - ((1-n)/(1+n))^2)
We used zeroang instead in Matlab which is probably why its in pvlib-python

@@ -831,7 +834,7 @@ def physicaliam(aoi, n=1.526, K=4., L=0.002):

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use
tau0 = (np.exp(- 1.0 * (K*L)) * (1 - (1 - n)^2/(1 + n)^2)
instead of the expression on line 826. I'd switch to tau0 when thetar_deg<zeroing.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with your suggestions, especially after having consulted my D&B book. I would suggest calculating reflectance and absorption separately before combining them in order to clarify the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I applied both your suggestions correctly. Could you take a look?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It wasn't quite correct, so I thought it might be easier to communicate it this way:

    # angle of refraction
    thetar_deg = tools.asind(tools.sind(aoi) / n)

    # reflectance and transmittance for normal incidence light
    rho_zero = ((1-n) / (1+n)) ** 2
    tau_zero = np.exp(-K*L)

    # reflectance for parallel and perpendicular polarized light
    rho_para = (tools.tand(thetar_deg - aoi) /
                tools.tand(thetar_deg + aoi)) ** 2
    rho_perp = (tools.sind(thetar_deg - aoi) /
                tools.sind(thetar_deg + aoi)) ** 2

    # transmittance for non-normal light
    tau = np.exp(-K*L / tools.cosd(thetar_deg))

    # iam is ratio of non-normal to normal incidence transmitted light
    # after deducting the reflected portion of each
    iam = ((1 - (rho_para + rho_perp) / 2) / (1 - rho_zero)
          * tau / tau_zero )


    # angles near zero produce nan, but iam is defined as one
    small_angle = 1e-06
    iam = np.where(np.abs(aoi) < small_angle, 1.0, iam)

    # angles at 90 degrees can produce tiny negative values, which should be zero
    # this is a result of calculation precision rather than the physical model
    iam = np.where(iam <  0, 0, iam)    

    # for light coming from behind the plane, none can enter the module
    iam = np.where(aoi > 90, 0, iam)

    # negative angles are not logical
    iam = np.where(aoi <  0, np.nan, iam)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the detailed explanation! So does this refer to reference 2 in the Docstring (I don't have access so cannot check)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ref [1] simply describes what is in ref [2].

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I implemented your suggestion except for setting negative angles to np.nan, to keep consistency with the other IAM functions for now...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consistency has some appeal...

@adriesse
Copy link
Member

adriesse commented Jul 6, 2017

Another observation: contrary to the statement in ref [1] the default value of 2 mm for the glass thickness is definitely at the thin extreme. The value I currently see most often is 3.2 mm. Previously I believe 4 mm was more common.

@jkfm
Copy link
Contributor Author

jkfm commented Jul 18, 2017

I agree with the default value, but I could not find any sources on it. Do you think, we should change the default nevertheless?

@adriesse
Copy link
Member

The source for the 2mm default has no sources either...

@cwhanse
Copy link
Member

cwhanse commented Jul 18, 2017

If it helps, the model output is essentially insensitive to the thickness parameter. You can double from 0.002 to 0.004 and the IAM changes very slightly.

For myself, I'd prefer the default to remain the value in the reference paper, simply for traceability back to literature. I do see merit in changing to 3.2mm, though.

@wholmgren
Copy link
Member

I suggest keeping the default set to 2 mm. I'll merge tomorrow unless there are objections.

@wholmgren wholmgren merged commit 3b0a86b into pvlib:master Aug 1, 2017
@jkfm jkfm deleted the modelchain_zero_instead_of_npnan branch December 14, 2017 13:25
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants