Changeset 282

Timestamp:

09/12/07 17:43:11 (17 years ago)

Author:

smasson

Message:

use of double precision to avoid rounding errors

Location:

trunk/SRC/Interpolation

Files:

• compute_fromirr_bilinear_weigaddr.pro (modified) (3 diffs)
• compute_fromreg_bilinear_weigaddr.pro (modified) (1 diff)
• compute_fromreg_imoms3_weigaddr.pro (modified) (4 diffs)
• quadrilateral2square.pro (modified) (3 diffs)
• square2quadrilateral.pro (modified) (2 diffs)

trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

@@ -216 +216 @@
                                     , newcoord[0, 2], newcoord[1, 2] $
                                     , newcoord[0, 1], newcoord[1, 1] $
-                                    , newcoord[0, 4], newcoord[1, 4])
+                                    , newcoord[0, 4], newcoord[1, 4], /double)
         ENDIF ELSE BEGIN
           xy = quadrilateral2square(xcell[0, ind], ycell[0, ind] $
                                     , xcell[1, ind], ycell[1, ind] $
                                     , xcell[2, ind], ycell[2, ind] $
-                                    , xcell[3, ind], ycell[3, ind], xx, yy)
+                                    , xcell[3, ind], ycell[3, ind], xx, yy, /double)
         ENDELSE
 ; take care of rounding errors...
-        zero = where(abs(xy) LT 1e-4)
-        IF zero[0] NE -1 THEN xy[zero] = 0
-        one = where(abs(1-xy) LT 1e-4)
-        IF one[0] NE -1 THEN xy[one] = 1
+        zero = where(abs(xy) LT 1e-6)
+        IF zero[0] NE -1 THEN xy[zero] = 0.d
+        one = where(abs(1.d - xy) LT 1e-6)
+        IF one[0] NE -1 THEN xy[one] = 1.d
 ; some checks...
         tmpmsk = omsk[celladdr[*, ind]]
@@ -249 +249 @@
           END
         ENDCASE
-
-
-
       ENDIF ELSE foraddr[n] = -1
     ENDIF ELSE foraddr[n] = -1
@@ -289 +286 @@
 ; mask the weight to suppress the corner located on land
   newaweig = newaweig*((omsk)[newaaddr])
-  totalweig = total(newaweig, 1)
 ; for cell with some land corner,
-; we have to redistribute the weight on the reaining water corners
+; we have to redistribute the weight on the remaining water corners
 ; weights normalization
-  totalweig = total(newaweig, 1)
-  IF min(totalweig, max = ma) EQ 0 then stop
-  IF ma GT 1 then stop
-  newaweig = newaweig/(replicate(1., 4)#totalweig)
-  totalweig = total(newaweig, 1)
-
+  totalweig = total(newaweig, 1, /double)
+  IF abs(min(totalweig, max = ma)) LT 1.e-6 then stop
+  IF abs(1.d - ma) GT 1.e-6 then stop
+  newaweig = newaweig/(replicate(1.d, 4)#totalweig)
 ; weights
   weig = dblarr(4, jpia*jpja)

trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro

@@ -217 +217 @@
   ENDIF
 ; check totalweight = 1
-  totalweig = abs(1-total(weig, 1))
+  totalweig = abs(1.d - total(weig, 1, /double))
   IF (where(temporary(totalweig) GE 1.e-5))[0] NE -1 THEN stop
 ;

trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro

@@ -284 +284 @@
   weig[15, short] = wx3 * wy3
 ;
-  mi = min(total(weig[*, short], 1), max = ma)
+  mi = min(total(weig[*, short], 1, /double), max = ma)
   IF abs(mi-1) GE 1.e-6 THEN stop
   IF abs(ma-1) GE 1.e-6 THEN stop
@@ -385 +385 @@
     weig[7, short] = wx3 * wy2
 ;
-    mi = min(total(weig[*, short], 1), max = ma)
+    mi = min(total(weig[*, short], 1, /double), max = ma)
     IF abs(mi-1) GE 1.e-6 THEN stop
     IF abs(ma-1) GE 1.e-6 THEN stop
@@ -451 +451 @@
     weig[3, short] = wx3
 ;
-    mi = min(total(weig[*, short], 1), max = ma)
+    mi = min(total(weig[*, short], 1, /double), max = ma)
     IF abs(mi-1) GE 1.e-6 THEN stop
     IF abs(ma-1) GE 1.e-6 THEN stop
@@ -517 +517 @@
     weig[3, short] = wx3
 ;
-    mi = min(total(weig[*, short], 1), max = ma)
+    mi = min(total(weig[*, short], 1, /double), max = ma)
     IF abs(mi-1) GE 1.e-6 THEN stop
     IF abs(ma-1) GE 1.e-6 THEN stop

trunk/SRC/Interpolation/quadrilateral2square.pro

@@ -38 +38 @@
 ; Can be scalar or array.
 ;
-; @keyword PERF
+; @keyword PERF {type=salar 0 or 1}{default=0}
+; activate to print the elapsed time spent within quadrilateral2square
+;
+; @keyword DOUBLE {type=salar 0 or 1}{default=0}
+; activate to perform double precision computation
 ;
 ; @returns
@@ -75 +79 @@
 ;-
 ;
-FUNCTION quadrilateral2square, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, PERF = perf
+FUNCTION quadrilateral2square, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, PERF = perf, DOUBLE = double
 ;
   compile_opt idl2, strictarrsubs
@@ -111 +115 @@
 ; get the matrix A
 ;
-  a = square2quadrilateral(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in)
+  a = square2quadrilateral(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, DOUBLE = double)
 ;
 ; compute the adjoint matrix

trunk/SRC/Interpolation/square2quadrilateral.pro

@@ -29 +29 @@
 ;
 ; @param xxin {in}{optional}
+; first coordinates of the point(s) for which we want to do the mapping.
 ; @param yyin {in}{optional}
-; the coordinates of the point(s) for which we want to do the mapping.
+; second coordinates of the point(s) for which we want to do the mapping.
+;
+; @keyword DOUBLE {type=salar 0 or 1}{default=0}
+; activate to perform double precision computation
 ;
 ; @returns
@@ -67 +71 @@
 ;-
 ;
-FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin
+FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, DOUBLE = double
 ;
 ; Warning, wrong definition of (x2,y2) and (x3,y3) at the bottom of