Earthdata Forum

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 2014-12-26 15:04
and later adapted to matlab by Reisinger (2015-01-12 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 2015-01-22 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=90-40, solz=90-70, relaz=175-180, ws=7)
That is, I have to enter 90-40, 90-70 and 175-180 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

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

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

Well, you could try:
M. Wang and S. Bailey, "Correction of sun glint contamination on the SeaWiFS ocean and atmosphere products," Appl. Opt.  40, 4790-4798 (2001).
and references therein, particularly
H. Gordon and M. Wang, "Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS: a preliminary algorithm," Appl. Opt.  33, 443-452 (1994).
which references
H. Gordon and M. Wang, "Surface-roughness considerations for atmospheric correction of ocean color sensors. 1: The Rayleigh-scattering component," Appl. Opt.  31, 4247-4260 (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

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?

Looks like the calculation of alpha in getglint.f is unnecessary with the assumptions made.

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.

I assume the code (line 59) is correct, ...

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 so

Sean

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