Updated working copy to latest revision, merged changes

.gitignore

@@ -35,3 +35,4 @@ nosetests.xml
 .idea
 *.iml
 *.ipr
+*~

gaip/Makefile

@@ -1,12 +1,12 @@
-FC := gfortran
+FC := ifort
 F2PY:= f2py
 
 all: _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so unmiximage.pyf unmiximage.so _interpolation.so
 
 _cast_shadow_mask.so _slope_self_shadow.so _surface_reflectance.so sys_variables.o sys_variables.mod _sat_sol_angles.so _satellite_model.so _track_time_info.so:
 	$$(module load gcc/5.2.0)
 	$(FC) -c sys_variables.f90
-	python setup.py build_ext -i
+	python setup.py config_fc --opt "-O3 -g" --arch "-xHost" build_ext -i
 
 unmiximage.pyf unmiximage.so:
 	$(F2PY) -h unmiximage.pyf -m unmiximage unmiximage.f90

gaip/black_sky.f90

@@ -1,6 +1,8 @@
-real function black_sky(brdf0,brdf1,brdf2,theta)
+function black_sky(brdf0,brdf1,brdf2,theta,veclen)
+  integer veclen
       real brdf0,brdf1,brdf2
-      real theta
+      real theta(veclen)
+      real black_sky(veclen)
 
 !   theta should be in radians
     black_sky = brdf0 + &

gaip/brdf_shape.f90

@@ -1,51 +1,71 @@
-real function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2)
+function RL_brdf (solar,view,ra,hb,br,brdf0,brdf1,brdf2,veclen)
 !   calculate the normalised BRDF shape function
-    real pi
-    real solar, view, ra, hb, br
+    integer veclen
+    real solar(veclen), view(veclen), ra(veclen), hb, br
     real brdf0,brdf1,brdf2
-    real rs_thick,li_sparse,secsolar,secvia
-    real cossolar,cosvia,cosra,sinsolar,sinvia,sinra
-    real cosxi,xi,tansolar,tanvia,theta_new_v,theta_new_s
-    real d_li2,x_li,cosl,l_li,o_li
+    real RL_brdf(veclen)
+
+    real rs_thick(veclen),li_sparse(veclen),secsolar(veclen),secvia(veclen)
+    real cossolar(veclen),cosvia(veclen),cosra(veclen),sinsolar(veclen),sinvia(veclen),sinra(veclen)
+    real cosxi(veclen),xi(veclen),tansolar(veclen),tanvia(veclen),theta_new_v(veclen),theta_new_s(veclen)
+    real d_li2(veclen),x_li(veclen),cosl(veclen),l_li(veclen),o_li(veclen)
+
+    real tan_theta_new_v(veclen), cos_theta_new_v(veclen), sin_theta_new_v(veclen)
+    real tan_theta_new_s(veclen), cos_theta_new_s(veclen), sin_theta_new_s(veclen)
+
+    real, parameter :: pi = 3.14159265358979
 
     RL_brdf=0.0
 
-    pi = 4 * atan(1.0)
+!    pi = 4 * atan(1.0)
 !   calculate Ross-thick kernel
     cossolar = cos(solar)
     cosvia = cos(view)
     cosra = cos(ra)
     sinsolar = sin(solar)
     sinvia = sin(view)
     sinra = sin(ra)
-    cosxi = cossolar*cosvia+sinsolar*sinvia*cosra
-    if (cosxi .ge.1) then
-        cosxi = 1
-    endif
+    cosxi = min(cossolar*cosvia+sinsolar*sinvia*cosra,1.0)
+!    if (cosxi .ge.1) then
+!        cosxi = 1
+!    endif
     xi = acos(cosxi)
-    rs_thick = ((pi/2.0-xi)*cos(xi)+sin(xi)) / (cossolar+cosvia) - pi/4.0
+    rs_thick = ((pi/2.0-xi)*cosxi+sin(xi)) / (cossolar+cosvia) - pi/4.0
 
 !   alculate Li-sparse
     tansolar = sinsolar/cossolar
     tanvia = sinvia/cosvia
-    theta_new_v = atan(br*tanvia)
-    theta_new_s = atan(br*tansolar)
-    cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra
-    if (cosxi.ge.1) then
-        cosxi = 1
-    endif
-    secsolar = 1.0/cos(theta_new_s)
-    secvia = 1.0/cos(theta_new_v)
-    d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra)
-    x_li = tan(theta_new_s)*tan(theta_new_v)*sinra
+!    theta_new_v = atan(br*tanvia)
+!    theta_new_s = atan(br*tansolar)
+    tan_theta_new_v = br*tanvia
+    ! Trig identities
+    cos_theta_new_v = 1.0 / sqrt(1 + tan_theta_new_v * tan_theta_new_v)
+    sin_theta_new_v = tan_theta_new_v * cos_theta_new_v
+
+    tan_theta_new_s = br*tansolar
+    ! Trig identities
+    cos_theta_new_s = 1.0 / sqrt(1 + tan_theta_new_s * tan_theta_new_s)
+    sin_theta_new_s = tan_theta_new_s * cos_theta_new_s
+
+!    cosxi = cos(theta_new_s)*cos(theta_new_v)+sin(theta_new_s)*sin(theta_new_v)*cosra
+    cosxi = min( cos_theta_new_s*cos_theta_new_v + sin_theta_new_s*sin_theta_new_s*cosra, 1.0 )
+!    if (cosxi.ge.1) then
+!        cosxi = 1
+!    endif
+    secsolar = 1.0/cos_theta_new_s
+    secvia = 1.0/cos_theta_new_v
+!    d_li2 = abs(tan(theta_new_s)**2+tan(theta_new_v)**2 - 2.0*tan(theta_new_s)*tan(theta_new_v)*cosra)
+    d_li2 = abs(tan_theta_new_s**2+tan_theta_new_v**2 - 2.0*tan_theta_new_s*tan_theta_new_v*cosra)
+    x_li = tan_theta_new_s*tan_theta_new_v*sinra
     cosl = hb*sqrt(d_li2+x_li**2)/(secsolar+secvia)
-    if (cosl .ge. 1.0) then
-        o_li=0.0
-    else
-        l_li=acos(cosl)
-        o_li=(l_li-sin(l_li)*cos(l_li))*(secsolar+secvia)/pi
-    endif
-       li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia
-        RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse
+!    if (cosl .ge. 1.0) then
+!    o_li=0.0
+!    else
+    l_li=acos(cosl)
+    o_li=(l_li-sin(l_li)*cosl)*(secsolar+secvia)/pi
+    where( cosl .ge. 1.0 ) o_li = 0.0
+!    endif
+    li_sparse=o_li-(secsolar+secvia)+0.5*(1.0+cosxi) * secsolar*secvia
+    RL_brdf=brdf0+brdf1*rs_thick+brdf2*li_sparse
     return
 end

gaip/setup.py

@@ -25,9 +25,7 @@
                                       'geo2metres_pixel_size.f90']),
                    Extension(name='_surface_reflectance',
                              sources=['surface_reflectance.f90',
-                                      'white_sky.f90',
-                                      'black_sky.f90',
-                                      'brdf_shape.f90']),
+                                      'reflectance_mod.f90']),
                    Extension(name='_satellite_model',
                              sources=['geo2metres_pixel_size.f90',
                                       'satellite_model.f90']),