Changeset 247 for codes

Timestamp:

07/22/14 16:27:58 (10 years ago)

Author:

dubos

Message:

Numerically stable centroid/circumcenter

File:

• codes/icosagcm/trunk/src/spherical_geom.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/spherical_geom.f90

@@ -179 +179 @@
 
 
-  SUBROUTINE circumcenter(A0,B0,C0,Center)
+  SUBROUTINE circumcenter(A0,B0,C0,center)
   USE vector
   IMPLICIT NONE
@@ -185 +185 @@
     REAL(rstd), INTENT(OUT) :: Center(3)
     
-    REAL(rstd)  :: a(3),b(3),c(3)
+    REAL(rstd)  :: a(3),b(3),c(3), ac(3), ab(3), p1(3), q(3), p2(3)
     
     a=A0/sqrt(sum(A0**2))
@@ -191 +191 @@
     c=C0/sqrt(sum(C0**2))
     
-    CALL Cross_product2(b-a,c-b,center)
-    center=center/sqrt(sum(center**2))
-    
+    ab=b-a
+    ac=c-a
+    CALL Cross_product2(ab,ac,p1)
+    IF(.FALSE.) THEN ! Direct solution, round-off error
+       center=p1/norm(p1)
+    ELSE ! Two-step solution, stable
+       q = SUM(ac**2)*ab-SUM(ab**2)*ac
+       CALL Cross_product2(p1,q,p2)
+       p2 = a + p2/(2.*SUM(p1**2))
+       center = p2/norm(p2)
+    END IF
   END SUBROUTINE circumcenter
     
@@ -204 +212 @@
     REAL(rstd), INTENT(OUT) :: Centr(3)
     
-    REAL(rstd) :: p1(3),p2(3),cross(3)
-    REAL(rstd) :: norm_cross
+    REAL(rstd) :: p1(3),p2(3),cross(3), cc(3)
+    REAL(rstd) :: norm_cross, area
     INTEGER :: i,j
-      
-      Centr(:)=0
-      DO i=1,n
-        j=MOD(i,n)+1
-        p1=points(:,i)/norm(points(:,i))
-        p2=points(:,j)/norm(points(:,j))
-        CALL cross_product2(p1,p2,cross)
-        norm_cross=norm(cross)
-        if (norm_cross<1e-10) CYCLE 
-        
-        Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross
-      ENDDO
-      
-      Centr(:)=centr(:)/norm(centr(:))
-  
+
+    centr(:)=0
+    IF(.FALSE.) THEN 
+       ! Gauss formula (subject to round-off error)
+       DO i=1,n
+          j=MOD(i,n)+1
+          p1=points(:,i)/norm(points(:,i))
+          p2=points(:,j)/norm(points(:,j))
+          CALL cross_product2(p1,p2,cross)
+          norm_cross=norm(cross)
+          if (norm_cross<1e-10) CYCLE 
+          centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross
+       ENDDO       
+    ELSE
+       ! Simple area-weighted average (second-order accurate, stable)
+       cc=SUM(points,2) ! arithmetic average used as center
+       cc=cc/norm(cc)
+       DO i=1,n
+          j=MOD(i,n)+1
+          p1=points(:,i)/norm(points(:,i))
+          p2=points(:,j)/norm(points(:,j))
+          CALL surf_triangle(cc,p1,p2,area)
+          centr(:)=centr(:)+area*(p1+p2+cc)
+       ENDDO
+    END IF
+    
+    centr(:)=centr(:)/norm(centr(:))
+
   END SUBROUTINE compute_centroid