Sunglint calculation
Sunglint calculation
ref. https://oceancolor.gsfc.nasa.gov/forum/oceancolor/topic_show.pl?tid=5643
(which is locked)
I apologize for raising this old thread, but I would like to predict sunglint using the code sub_glint.pl contributed by Sean Bailey on 20141226 15:04
and later adapted to matlab by Reisinger (20150112 17:33 ).
I've reproduced this code (in R) but I'm a bit confused by the comments at the end of the thread and would like to clarify.
I understand I have to keep:
senz = deg2rad(senz);
solz = deg2rad(solz);
relaz = deg2rad(relaz);
(see message By SeanBailey Date 20150122 08:19 )
but in that case, I understand that, in glint_example.txt:
1. the values of senz,solz are not zenith angles, but elevation angles,
2. the relaz value has to be converted from i.e. 175 to 175 180
In that case, I get the same value of g as in glint_example.txt for 40, 70, 175
with the following input for my R code
glint(senz=9040, solz=9070, relaz=175180, ws=7)
That is, I have to enter 9040, 9070 and 175180 to get the same output as recorded in glint_example.txt for 40, 70 and 175.
To make sure I'm doing right, could you confirm the above 2 questions? What makes me doubt is that if senz and solz are,
in fact, elevation angles,
your example glint_example.txt calculates glint for sensor elevations of 0 to 45, which seems really weird.
3. I understand your code requires calculating relaz as
Relative Azimuth = Sensor Azimuth  180.0  Solar azimuth
as I actually do.
But I'm also a bit confused with the 180 in the relaz formula
(even after reading some other old discussions in this forum).
If a given pixel of seawifs has a sensor azimuth of 260, is the sensor looking west or east?
If the solar azimuth is 260, the sensor azimut for which the sensor
would look into the sun would be 260+180, correct? That way relaz will be 0 (= 260+180 180  260)
Finally
4. Is there a threshold of acceptabe g?
5. Any article behind this code for g?
Thanks
Agus
(which is locked)
I apologize for raising this old thread, but I would like to predict sunglint using the code sub_glint.pl contributed by Sean Bailey on 20141226 15:04
and later adapted to matlab by Reisinger (20150112 17:33 ).
I've reproduced this code (in R) but I'm a bit confused by the comments at the end of the thread and would like to clarify.
I understand I have to keep:
senz = deg2rad(senz);
solz = deg2rad(solz);
relaz = deg2rad(relaz);
(see message By SeanBailey Date 20150122 08:19 )
but in that case, I understand that, in glint_example.txt:
1. the values of senz,solz are not zenith angles, but elevation angles,
2. the relaz value has to be converted from i.e. 175 to 175 180
In that case, I get the same value of g as in glint_example.txt for 40, 70, 175
with the following input for my R code
glint(senz=9040, solz=9070, relaz=175180, ws=7)
That is, I have to enter 9040, 9070 and 175180 to get the same output as recorded in glint_example.txt for 40, 70 and 175.
To make sure I'm doing right, could you confirm the above 2 questions? What makes me doubt is that if senz and solz are,
in fact, elevation angles,
your example glint_example.txt calculates glint for sensor elevations of 0 to 45, which seems really weird.
3. I understand your code requires calculating relaz as
Relative Azimuth = Sensor Azimuth  180.0  Solar azimuth
as I actually do.
But I'm also a bit confused with the 180 in the relaz formula
(even after reading some other old discussions in this forum).
If a given pixel of seawifs has a sensor azimuth of 260, is the sensor looking west or east?
If the solar azimuth is 260, the sensor azimut for which the sensor
would look into the sun would be 260+180, correct? That way relaz will be 0 (= 260+180 180  260)
Finally
4. Is there a threshold of acceptabe g?
5. Any article behind this code for g?
Thanks
Agus
Filters:

 User Services
 Posts: 1483
 Joined: Wed Sep 18, 2019 6:15 pm America/New_York
 Been thanked: 5 times
Sunglint calculation
Agus,
The Perl code I originally posted was used with inputs coming from another program that output elevation, not zenith angles.
The 180 is to put the view in the direction of the glint not the direction of the sun. A relaz of 180 would be the sun behind the sensor and high glint.
We use a glint threshold of 0.001 for moderate glint and 0.005 for high glint.
Sean
The Perl code I originally posted was used with inputs coming from another program that output elevation, not zenith angles.
The 180 is to put the view in the direction of the glint not the direction of the sun. A relaz of 180 would be the sun behind the sensor and high glint.
We use a glint threshold of 0.001 for moderate glint and 0.005 for high glint.
Sean
Sunglint calculation
Thanks Sean,
any reference that inspired the code? Any way I can make reference
to your code or should I just put pers. comm. ?
Agus
any reference that inspired the code? Any way I can make reference
to your code or should I just put pers. comm. ?
Agus

 User Services
 Posts: 1483
 Joined: Wed Sep 18, 2019 6:15 pm America/New_York
 Been thanked: 5 times
Sunglint calculation
Well, you could try:
which references
which references
...but ultimately, I just transposed the FORTRAN code
Sean
M. Wang and S. Bailey, "Correction of sun glint contamination on the SeaWiFS ocean and atmosphere products," Appl. Opt. 40, 47904798 (2001).
and references therein, particularlyH. Gordon and M. Wang, "Retrieval of waterleaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm," Appl. Opt. 33, 443452 (1994).
which references
H. Gordon and M. Wang, "Surfaceroughness considerations for atmospheric correction of ocean color sensors. 1: The Rayleighscattering component," Appl. Opt. 31, 42474260 (1992).
which references
C. Cox and W. Munk, “Measurements of the roughness of the sea surface from photographs of the sun’s glitter,” J. Opt. Soc.Am. 44, 838 – 850 (1954).
...but ultimately, I just transposed the FORTRAN code
getglint.f
that is part of our atmospheric correction code into Perl...Sean

 Posts: 16
 Joined: Sun Jul 31, 2005 9:10 pm America/New_York
Re: Sunglint calculation
Looks like the calculation of alpha in getglint.f is unnecessary with the assumptions made.
When sigc = sigu,
expon = (swig**2+eta**2)/2
= ((sin(alphap)*tan(beta)/sigc)**2+(cos(alphap)*tan(beta)/sigu)**2)/2
= tan(beta)**2/sigc**2 * ((sin(alphap)**2+cos(alphap)**2)/2
= tan(beta)**2/sigc**2 /2
alphap got canceled out because the sum of its sine squared and cosine squared is always 1.
Then the calculation of slope probability distribution and normalized glint radiance can be made without any wind direction information, which is not needed as one of the inputs with the assumptions made in the simplification.
Relevant codes in getglint.f:
73 sigc = .04964*sqrt(y4)
74 sigu = .04964*sqrt(y4)
75 c endif
76 chi = x5
77 alphap = alpha+chi
78 swig = sin(alphap)*tan(beta)/sigc
79 eta = cos(alphap)*tan(beta)/sigu
80 expon = (swig**2+eta**2)/2.
81 if (expon.lt.30.) expon = 30. ! trap underflow
82 if (expon.gt.+30.) expon = +30. ! trap overflow
83 prob = exp(expon)/(2.*pi*sigu*sigc)
Another clarification question:
These two lines are inconsistent.
25 C cos(OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)
59 omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2.
I assume the code (line 59) is correct, and the correct formula on line 25 should be cos(2*OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)?
Eq. 3 in Cox and Munk (1954) is consistent with line 25 of getglint.f, also missing a factor of 2 in front of omega?
When sigc = sigu,
expon = (swig**2+eta**2)/2
= ((sin(alphap)*tan(beta)/sigc)**2+(cos(alphap)*tan(beta)/sigu)**2)/2
= tan(beta)**2/sigc**2 * ((sin(alphap)**2+cos(alphap)**2)/2
= tan(beta)**2/sigc**2 /2
alphap got canceled out because the sum of its sine squared and cosine squared is always 1.
Then the calculation of slope probability distribution and normalized glint radiance can be made without any wind direction information, which is not needed as one of the inputs with the assumptions made in the simplification.
Relevant codes in getglint.f:
73 sigc = .04964*sqrt(y4)
74 sigu = .04964*sqrt(y4)
75 c endif
76 chi = x5
77 alphap = alpha+chi
78 swig = sin(alphap)*tan(beta)/sigc
79 eta = cos(alphap)*tan(beta)/sigu
80 expon = (swig**2+eta**2)/2.
81 if (expon.lt.30.) expon = 30. ! trap underflow
82 if (expon.gt.+30.) expon = +30. ! trap overflow
83 prob = exp(expon)/(2.*pi*sigu*sigc)
Another clarification question:
These two lines are inconsistent.
25 C cos(OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)
59 omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2.
I assume the code (line 59) is correct, and the correct formula on line 25 should be cos(2*OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)?
Eq. 3 in Cox and Munk (1954) is consistent with line 25 of getglint.f, also missing a factor of 2 in front of omega?

 User Services
 Posts: 1483
 Joined: Wed Sep 18, 2019 6:15 pm America/New_York
 Been thanked: 5 times
Re: Sunglint calculation
Yes. The code is capable of supporting wind direction, but it is not used to do so. The windspeed is also assumed to be isotropic (because of the lack of wind direction information), so there are indeed some unnecessary calculations being made for the way the code is used.Looks like the calculation of alpha in getglint.f is unnecessary with the assumptions made.
I *think* the code is correct, but haven't had time to verify it. If it is NOT correct, it's been wrong for over 30 years ...since you raise the possibility, I'll find the time (or delegate someone else) to do soI assume the code (line 59) is correct, ...
Sean

 Posts: 1
 Joined: Wed Aug 09, 2023 11:38 am America/New_York
Re: Sunglint calculation
starbright wrote:
> Another clarification question:
> These two lines are inconsistent.
> 25 C cos(OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)
> 59 omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2.
>
> I assume the code (line 59) is correct, and the correct formula on line 25
> should be cos(2*OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)?
>
> Eq. 3 in Cox and Munk (1954) is consistent with line 25 of getglint.f, also
> missing a factor of 2 in front of omega?
The equations are correct.
The equation on line 59
omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2
is in fact Equation 8:
cos(2*OMEGA) = sin(PHI)sin(MU)cos(NU)+cos(PHI)cos(MU)
MU : y1, sensor zenith angle;
PHI : y2, solar zenith angle;
NU : y3 + pi, relative azimuth angle. Units are radians. Keep in mind that we define relative azimuth (y3) as
y3 = sensor_azimuth  pi  solar_azimuth.
and NU is "conventional" relative azimuth angle, NU = sensor_azimuth  solar_azimuth
thus
cos(NU) =  cos(y3)
2*OMEGA is an angle between the incident ray and the reflected rays. It can be computed from a dot product between these two vectors which gives Eq. 8 showed above.
Equation 3 in the Cox and Munk paper indeed corresponds to the equation on line 25 in the code but it's uses BETA, which, in fact, computed later in the same code (line 61):
beta = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
Thanks
Alex
> Another clarification question:
> These two lines are inconsistent.
> 25 C cos(OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)
> 59 omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2.
>
> I assume the code (line 59) is correct, and the correct formula on line 25
> should be cos(2*OMEGA) = cos(BETA)cos(PHI)sin(BETA)sin(PHI)cos(ALPHA)?
>
> Eq. 3 in Cox and Munk (1954) is consistent with line 25 of getglint.f, also
> missing a factor of 2 in front of omega?
The equations are correct.
The equation on line 59
omega = acoss(cos(y1)*cos(y2)sin(y1)*sin(y2)*cos(y3))/2
is in fact Equation 8:
cos(2*OMEGA) = sin(PHI)sin(MU)cos(NU)+cos(PHI)cos(MU)
MU : y1, sensor zenith angle;
PHI : y2, solar zenith angle;
NU : y3 + pi, relative azimuth angle. Units are radians. Keep in mind that we define relative azimuth (y3) as
y3 = sensor_azimuth  pi  solar_azimuth.
and NU is "conventional" relative azimuth angle, NU = sensor_azimuth  solar_azimuth
thus
cos(NU) =  cos(y3)
2*OMEGA is an angle between the incident ray and the reflected rays. It can be computed from a dot product between these two vectors which gives Eq. 8 showed above.
Equation 3 in the Cox and Munk paper indeed corresponds to the equation on line 25 in the code but it's uses BETA, which, in fact, computed later in the same code (line 61):
beta = acoss((cos(y1)+cos(y2))/(2.*cos(omega)))
Thanks
Alex