Okay, hot on the heels of yesterday's "how come ccdpack has forked kappa
routines" discussions here comes the rematch of "how come kappa has loads
of ccdpack routines".
KAPPA has 75 ccdpack routines, all but one are derived from generic source
without actually having the generic source in CVS [ccdpack does]
The diffs are below, the main difference seems to be that Dave has added
extra IMETH switches to support for unweighted mean and also support for
IMETH=300 and an alternative way of calculating variance by estimation.
ccg1_mdr3 seems to be identical except for a missing bug fix comment from
Mark (it seems that David fixed it during the same period).
David's additions seem to be plausible for CCDPACK (a couple would
require argument tweaks but nothing that couldn't be specified as a
constant).
--
Tim Jenness
JAC software
http://www.jach.hawaii.edu/~timj
Fortran file: ccd1_orva Generic: NONE
--- ccg1_cm1d.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_cm1dd.f 2004-09-07 10:03:24.000000000 -1000
@@ -41,7 +41,8 @@
* IMETH = INTEGER (Given)
* The method to use in combining the lines. Has a code of 2 to 10
* which represent.
-* 2 = MEAN
+* 1 = UNWEIGHTED MEAN
+* 2 = WEIGHTED MEAN
* 3 = MEDIAN
* 4 = TRIMMED MEAN
* 5 = MODE
@@ -50,6 +51,7 @@
* 8 = MINMAX MEAN
* 9 = BROADENED MEDIAN
* 10 = SIGMA CLIPPED MEDIAN
+* 300 = Median, but estimating variance from mean variance.
* MINPIX = INTEGER (Given)
* The minimum number of pixels required to contribute to an
* output pixel.
@@ -75,7 +77,7 @@
* Workspace for order statistics calculations.
* COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given and Returned)
* Workspace for storing ordered statistics variace-covariance
-* matrix.
+* matrix. Not used for IMETHs 1, 2 or 300.
* NCON( NLINES ) = DOUBLE PRECISION (Given and Returned)
* The actual number of contributing pixels from each input line
* to the output line.
@@ -123,6 +125,14 @@
* Added MEDCLIP method.
* 16-NOV-1998 (PDRAPER):
* Added fast median.
+* 14-DEC-2001 (DSB):
+* Added unweighted mean method.
+* 1-NOV-2002 (DSB):
+* Added workspace argument to CCD1_ORVAR.
+* 27-AUG-2003 (DSB):
+* Added IMETH 300 to provide a means of estimating the output
+* variances for large samples, for which the memory requirements
+* for the COVEC array would be too large.
* {enter_further_changes_here}
* Bugs:
@@ -135,6 +145,7 @@
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
+ INCLUDE 'CNF_PAR' ! For CNF_PVAL function
* Arguments Given:
INTEGER NPIX
@@ -166,6 +177,8 @@
* Status:
INTEGER STATUS ! Global status
+* Local Variables:
+ INTEGER IPW1 ! Work space pointer
*.
* Check inherited global status.
@@ -175,23 +188,33 @@
* for the order statistics of a normal popl with up to NLINE members.
* This also sets up the scale factor for converting mean variances to
* median variances.
- IF( IMETH .NE. 2 ) THEN
- CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC, STATUS )
+ IF( IMETH .NE. 1 .AND. IMETH .NE. 2 .AND. IMETH .NE. 300 ) THEN
+ CALL PSX_CALLOC( NLINES*NLINES, '_DOUBLE', IPW1, STATUS )
+ CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC,
+ : %VAL( CNF_PVAL( IPW1 ) ),
+ : STATUS )
+ CALL PSX_FREE( IPW1, STATUS )
END IF
* Now branch for each method.
- IF ( IMETH .EQ. 2 ) THEN
+ IF ( IMETH .EQ. 1 ) THEN
+
+* Forming the unweighted mean.
+ CALL CCG1_UMD1D( STACK, NPIX, NLINES, VARS, MINPIX,
+ : RESULT, RESVAR, NCON, STATUS )
+
+ ELSE IF ( IMETH .EQ. 2 ) THEN
* Forming the weighted mean.
CALL CCG1_MED1D( STACK, NPIX, NLINES, VARS, MINPIX,
: RESULT, RESVAR, NCON, STATUS )
- ELSE IF ( IMETH .EQ. 3 ) THEN
+ ELSE IF ( IMETH .EQ. 3 .OR. IMETH .EQ. 300 ) THEN
* Forming the weighted median.
- CALL CCG1_MDD1D( STACK, NPIX, NLINES, VARS, MINPIX, COVEC,
- : NMAT, RESULT, RESVAR, WRK1, WRK2, NCON,
- : POINT, USED, STATUS )
+ CALL CCG1_MDD1D( ( IMETH .EQ. 3 ), STACK, NPIX, NLINES, VARS,
+ : MINPIX, COVEC, NMAT, RESULT, RESVAR, WRK1,
+ : WRK2, NCON, POINT, USED, STATUS )
ELSE IF ( IMETH .EQ. 4 ) THEN
Fortran file: ccg1_cm1r Generic: ccg1_cm1r.gen
--- ccg1_cm1r.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_cm1rr.f 2004-09-07 10:03:24.000000000 -1000
@@ -41,7 +41,8 @@
* IMETH = INTEGER (Given)
* The method to use in combining the lines. Has a code of 2 to 9
* which represent.
-* 2 = MEAN
+* 1 = MEAN
+* 2 = WEIGHTED MEAN
* 3 = MEDIAN
* 4 = TRIMMED MEAN
* 5 = MODE
@@ -50,6 +51,7 @@
* 8 = MINMAX MEAN
* 9 = BROADENED MEDIAN
* 10 = SIGMA CLIPPED MEDIAN
+* 300 = Median, but estimating variance from mean variance.
* MINPIX = INTEGER (Given)
* The minimum number of pixels required to contribute to an
* output pixel.
@@ -75,7 +77,7 @@
* Workspace for order statistics calculations.
* COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given and Returned)
* Workspace for storing ordered statistics variace-covariance
-* matrix.
+* matrix. Not used for IMETHs 1, 2 or 300.
* NCON( NLINES ) = DOUBLE PRECISION (Given and Returned)
* The actual number of contributing pixels from each input line
* to the output line.
@@ -123,6 +125,14 @@
* Added sigma clipped median.
* 16-NOV-1998 (PDRAPER):
* Added fast median.
+* 9-SEP-2002 (DSB):
+* Added unweighted mean method.
+* 1-NOV-2002 (DSB):
+* Added workspace argument to CCD1_ORVAR.
+* 27-AUG-2003 (DSB):
+* Added IMETH 300 to provide a means of estimating the output
+* variances for large samples, for which the memory requirements
+* for the COVEC array would be too large.
* {enter_further_changes_here}
* Bugs:
@@ -135,6 +145,7 @@
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
+ INCLUDE 'CNF_PAR' ! For CNF_PVAL function
* Arguments Given:
INTEGER NPIX
@@ -166,6 +177,8 @@
* Status:
INTEGER STATUS ! Global status
+* Local Variables:
+ INTEGER IPW1 ! Work space pointer
*.
* Check inherited global status.
@@ -175,23 +188,33 @@
* for the order statistics of a normal popl with up to NLINE members.
* This also sets up the scale factor for converting mean variances to
* median variances.
- IF( IMETH .NE. 2 ) THEN
- CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC, STATUS )
+ IF( IMETH .NE. 1 .AND. IMETH .NE. 2 .AND. IMETH .NE. 300 ) THEN
+ CALL PSX_CALLOC( NLINES*NLINES, '_DOUBLE', IPW1, STATUS )
+ CALL CCD1_ORVAR( NLINES, NMAT, PP, COVEC,
+ : %VAL( CNF_PVAL( IPW1 ) ),
+ : STATUS )
+ CALL PSX_FREE( IPW1, STATUS )
END IF
* Now branch for each method.
- IF ( IMETH .EQ. 2 ) THEN
+ IF ( IMETH .EQ. 1 ) THEN
+
+* Forming the unweighted mean.
+ CALL CCG1_UMR1R( STACK, NPIX, NLINES, VARS, MINPIX,
+ : RESULT, RESVAR, NCON, STATUS )
+
+ ELSE IF ( IMETH .EQ. 2 ) THEN
* Forming the weighted mean.
CALL CCG1_MER1R( STACK, NPIX, NLINES, VARS, MINPIX,
: RESULT, RESVAR, NCON, STATUS )
- ELSE IF ( IMETH .EQ. 3 ) THEN
+ ELSE IF ( IMETH .EQ. 3 .OR. IMETH .EQ. 300 ) THEN
* Forming the weighted median.
- CALL CCG1_MDR1R( STACK, NPIX, NLINES, VARS, MINPIX, COVEC,
- : NMAT, RESULT, RESVAR, WRK1, WRK2, NCON,
- : POINT, USED, STATUS )
+ CALL CCG1_MDR1R( ( IMETH .EQ. 3 ), STACK, NPIX, NLINES, VARS,
+ : MINPIX, COVEC, NMAT, RESULT, RESVAR, WRK1,
+ : WRK2, NCON, POINT, USED, STATUS )
ELSE IF ( IMETH .EQ. 4 ) THEN
Fortran file: ccg1_cm3d Generic: ccg1_cm3d.gen
--- ccg1_cm3d.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_cm3dd.f 2002-09-10 00:12:55.000000000 -1000
@@ -39,7 +39,8 @@
* IMETH = INTEGER (Given)
* The method to use in combining the lines. Has a code of 2 to 10.
* which represent.
-* 2 = MEAN
+* 1 = MEAN
+* 2 = WEIGHTED MEAN
* 3 = MEDIAN
* 4 = TRIMMED MEAN
* 5 = MODE
@@ -100,6 +101,8 @@
* Added sigma clipped median.
* 16-NOV-1998 (PDRAPER):
* Added fast median.
+* 9-SEP-2002 (DSB):
+* Added unweighted mean method.
* {enter_further_changes_here}
* Bugs:
@@ -145,7 +148,13 @@
IF ( STATUS .NE. SAI__OK ) RETURN
* Branch for each method.
- IF ( IMETH .EQ. 2 ) THEN
+ IF ( IMETH .EQ. 1 ) THEN
+
+* Forming the unweighted mean.
+ CALL CCG1_UMD3D( STACK, NPIX, NLINES, VARS, MINPIX,
+ : RESULT, NCON, STATUS )
+
+ ELSE IF ( IMETH .EQ. 2 ) THEN
* Forming the weighted mean.
CALL CCG1_MED3D( STACK, NPIX, NLINES, VARS, MINPIX,
@@ -220,4 +229,4 @@
END IF
END
-* $Id: ccg1_cm3d.gen,v 1.3 1998/11/18 19:12:48 pdraper Exp $
+* $Id: ccg1_cm3dd.f,v 1.2 2002/09/10 10:12:55 dsb Exp $
Fortran file: ccg1_cm3r Generic: ccg1_cm3r.gen
--- ccg1_cm3r.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_cm3rr.f 2002-09-10 00:12:55.000000000 -1000
@@ -39,7 +39,8 @@
* IMETH = INTEGER (Given)
* The method to use in combining the lines. Has a code of 2 to 10.
* which represent.
-* 2 = MEAN
+* 1 = MEAN
+* 2 = WEIGHTED MEAN
* 3 = MEDIAN
* 4 = TRIMMED MEAN
* 5 = MODE
@@ -101,6 +102,8 @@
* Added sigma clipped median.
* 16-NOV-1998 (PDRAPER):
* Added fast median.
+* 9-SEP-2002 (DSB):
+* Added unweighted mean method.
* {enter_further_changes_here}
* Bugs:
@@ -146,7 +149,13 @@
IF ( STATUS .NE. SAI__OK ) RETURN
* Branch for each method.
- IF ( IMETH .EQ. 2 ) THEN
+ IF ( IMETH .EQ. 1 ) THEN
+
+* Forming the unweighted mean.
+ CALL CCG1_UMR3R( STACK, NPIX, NLINES, VARS, MINPIX,
+ : RESULT, NCON, STATUS )
+
+ ELSE IF ( IMETH .EQ. 2 ) THEN
* Forming the weighted mean.
CALL CCG1_MER3R( STACK, NPIX, NLINES, VARS, MINPIX,
Fortran file: ccg1_mdd1 Generic: ccg1_mdd1.gen
--- ccg1_mdd1.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_mdd1d.f 2003-08-27 08:20:06.000000000 -1000
@@ -1,6 +1,6 @@
- SUBROUTINE CCG1_MDD1D( STACK, NPIX, NLINES, VARS, MINPIX,
- : COVEC, NMAT, RESULT, RESVAR, WRK1,
- : WRK2, NCON, POINT, USED, STATUS )
+ SUBROUTINE CCG1_MDD1D( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
+ : COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
+ : NCON, POINT, USED, STATUS )
*+
* Name:
* CCG1_MDD1
@@ -12,9 +12,9 @@
* Starlink Fortran 77
* Invocation:
-* CALL CCG1_MDD1( STACK, NPIX, NLINES, VARS, MINPIX, COVEC, NMAT,
-* RESULT, RESVAR, WRK1, WRK2, NCON, POINT, USED,
-* STATUS )
+* CALL CCG1_MDD1D( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
+* COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
+* NCON, POINT, USED, STATUS )
* Description:
* This routine accepts an array consisting a series of (vectorised)
@@ -24,6 +24,11 @@
* in the array RESULT.
* Arguments:
+* CALCMV = LOGICAL (Given)
+* If .FALSE. then the output variances are estimated by scaling
+* the variance on the weighted mean (rather than the weighted median)
+* by Pi/2. Otherwise, the output variances are calculated using the
+* COVEC array.
* STACK( NPIX, NLINES ) = DOUBLE PRECISION (Given)
* The array of lines which are to be combined into a single line.
* NPIX = INTEGER (Given)
@@ -38,7 +43,7 @@
* COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given)
* The packed variance-covariance matrix of the order statistics
* from a normal distribution of sizes up to NLINES, produced by
-* CCD1_ORVAR.
+* CCD1_ORVAR. Not used if CALCMV is .FALSE.
* NMAT = INTEGER (Given)
* Size of the first dimension of COVEC.
* RESULT( NPIX ) = DOUBLE PRECISION (Returned)
@@ -67,6 +72,8 @@
* History:
* 21-MAY-1992 (PDRAPER):
* Original version.
+* 27-AUG-2003 (DSB):
+* Added argument CALCMV.
* {enter_changes_here}
* Bugs:
@@ -82,6 +89,7 @@
INCLUDE 'PRM_PAR' ! PRIMDAT constants
* Arguments Given:
+ LOGICAL CALCMV
INTEGER NPIX
INTEGER NLINES
INTEGER MINPIX
@@ -178,8 +186,13 @@
CALL CCG1_IS3D( WRK1, WRK2, POINT, NGOOD, STATUS )
* Find the weighted median.
- CALL CCG1_WTM3D( WRK1, WRK2, SVAR, NGOOD, USED,
- : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
+ IF( CALCMV ) THEN
+ CALL CCG1_WTM3D( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
+ : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
+ ELSE
+ CALL CCG1_WTM3D( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
+ : COVEC, VAL, VAR, STATUS )
+ END IF
* Update the used line counters.
DO 3 J = 1, NGOOD
@@ -211,4 +224,4 @@
CALL NUM_REVRT
END
-* $Id: ccg1_mdd1.gen,v 1.3 1998/11/20 16:05:26 pdraper Exp $
+* $Id: ccg1_mdd1d.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
Fortran file: ccg1_mdd3 Generic: ccg1_mdd3.gen
--- ccg1_mdd3.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_mdd3d.f 2000-09-05 03:02:58.000000000 -1000
@@ -53,15 +53,11 @@
* Authors:
* PDRAPER: Peter Draper (STARLINK)
-* MBT: Mark Taylor (STARLINK)
* {enter_new_authors_here}
* History:
* 19-MAY-1992 (PDRAPER):
* Original version.
-* 27-NOV-2000 (MBT):
-* Fixed bug; variance values were not sorted corresponding to
-* order of data values being stacked.
* {enter_changes_here}
* Bugs:
@@ -157,7 +153,7 @@
IF ( NGOOD .GT. 0 ) THEN
* Sort these values into increasing order.
- CALL CCG1_IS3D( WRK1, WRK2, POINT, NGOOD, STATUS )
+ CALL CCG1_IS2D( WRK1, POINT, NGOOD, STATUS )
* Find the weighted median.
CALL CCG1_WTM2D( WRK1, WRK2, NGOOD, USED, VAL, STATUS )
@@ -189,4 +185,4 @@
CALL NUM_REVRT
END
-* $Id: ccg1_mdd3.gen,v 1.5 2000/11/28 13:12:14 mbt Exp $
+* $Id: ccg1_mdd3d.f,v 1.1 2000/09/05 13:02:58 dsb Exp $
Fortran file: ccg1_mdr1 Generic: ccg1_mdr1.gen
--- ccg1_mdr1.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_mdr1r.f 2003-08-27 08:20:06.000000000 -1000
@@ -1,6 +1,6 @@
- SUBROUTINE CCG1_MDR1R( STACK, NPIX, NLINES, VARS, MINPIX,
- : COVEC, NMAT, RESULT, RESVAR, WRK1,
- : WRK2, NCON, POINT, USED, STATUS )
+ SUBROUTINE CCG1_MDR1R( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
+ : COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
+ : NCON, POINT, USED, STATUS )
*+
* Name:
* CCG1_MDR1
@@ -12,9 +12,9 @@
* Starlink Fortran 77
* Invocation:
-* CALL CCG1_MDR1( STACK, NPIX, NLINES, VARS, MINPIX, COVEC, NMAT,
-* RESULT, RESVAR, WRK1, WRK2, NCON, POINT, USED,
-* STATUS )
+* CALL CCG1_MDR1R( CALCMV, STACK, NPIX, NLINES, VARS, MINPIX,
+* COVEC, NMAT, RESULT, RESVAR, WRK1, WRK2,
+* NCON, POINT, USED, STATUS )
* Description:
* This routine accepts an array consisting a series of (vectorised)
@@ -24,6 +24,11 @@
* output means are returned in the array RESULT.
* Arguments:
+* CALCMV = LOGICAL (Given)
+* If .FALSE. then the output variances are estimated by scaling
+* the variance on the weighted mean (rather than the weighted median)
+* by Pi/2. Otherwise, the output variances are calculated using the
+* COVEC array.
* STACK( NPIX, NLINES ) = REAL (Given)
* The array of lines which are to be combined into a single line.
* NPIX = INTEGER (Given)
@@ -38,7 +43,7 @@
* COVEC( NMAT, NLINES ) = DOUBLE PRECISION (Given)
* The packed variance-covariance matrix of the order statistics
* from a normal distribution of sizes up to NLINES, produced by
-* CCD1_ORVAR.
+* CCD1_ORVAR. Not used if CALCMV is .FALSE.
* NMAT = INTEGER (Given)
* Size of the first dimension of COVEC.
* RESULT( NPIX ) = REAL (Returned)
@@ -67,6 +72,8 @@
* History:
* 20-MAY-1992 (PDRAPER):
* Original version.
+* 27-AUG-2003 (DSB):
+* Added argument CALCMV.
* {enter_changes_here}
* Bugs:
@@ -82,6 +89,7 @@
INCLUDE 'PRM_PAR' ! PRIMDAT constants
* Arguments Given:
+ LOGICAL CALCMV
INTEGER NPIX
INTEGER NLINES
INTEGER MINPIX
@@ -178,8 +186,13 @@
CALL CCG1_IS3R( WRK1, WRK2, POINT, NGOOD, STATUS )
* Find the weighted median.
- CALL CCG1_WTM3R( WRK1, WRK2, SVAR, NGOOD, USED,
- : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
+ IF( CALCMV ) THEN
+ CALL CCG1_WTM3R( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
+ : COVEC( 1, NGOOD ), VAL, VAR, STATUS )
+ ELSE
+ CALL CCG1_WTM3R( CALCMV, WRK1, WRK2, SVAR, NGOOD, USED,
+ : COVEC, VAL, VAR, STATUS )
+ END IF
* Update the used line counters.
DO 3 J = 1, NGOOD
@@ -211,4 +224,4 @@
CALL NUM_REVRT
END
-* $Id: ccg1_mdr1.gen,v 1.3 1998/11/20 16:05:26 pdraper Exp $
+* $Id: ccg1_mdr1r.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
Fortran file: ccg1_mdr3 Generic: ccg1_mdr3.gen
--- ccg1_mdr3.f 2004-10-04 16:25:59.000000000 -1000
+++ ccg1_mdr3r.f 2001-11-13 03:32:18.000000000 -1000
@@ -53,15 +53,11 @@
* Authors:
* PDRAPER: Peter Draper (STARLINK)
-* MBT: Mark Taylor (STARLINK)
* {enter_new_authors_here}
* History:
* 19-MAY-1992 (PDRAPER):
* Original version.
-* 27-NOV-2000 (MBT):
-* Fixed bug; variance values were not sorted corresponding to
-* order of data values being stacked.
* {enter_changes_here}
* Bugs:
@@ -189,4 +185,4 @@
CALL NUM_REVRT
END
-* $Id: ccg1_mdr3.gen,v 1.5 2000/11/28 13:12:14 mbt Exp $
+* $Id: ccg1_mdr3r.f,v 1.2 2001/11/13 13:32:18 dsb Exp $
--- ccg1_wtm3.f 2004-10-04 16:26:01.000000000 -1000
+++ ccg1_wtm3d.f 2003-08-27 08:20:06.000000000 -1000
@@ -1,5 +1,5 @@
- SUBROUTINE CCG1_WTM3D( ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
- : RESULT, RESVAR, STATUS )
+ SUBROUTINE CCG1_WTM3D( CALCMV, ORDDAT, WEIGHT, VAR, NENT, USED,
+ : COVAR, RESULT, RESVAR, STATUS )
*+
* Name:
* CCG1_WTM3D
@@ -13,8 +13,8 @@
* Starlink Fortran 77
* Invocation:
-* CALL CCG1_WTM3D( ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
-* RESULT, RESVAR, STATUS )
+* CALL CCG1_WTM3D( CALCMV, ORDDAT, WEIGHT, VAR, NENT, USED, COVAR,
+* RESULT, RESVAR, STATUS )
* Description:
* This routine finds a value which can be associated with the half-
@@ -33,6 +33,11 @@
* variance.
* Arguments:
+* CALCMV = LOGICAL (Given)
+* If .FALSE. then the output variances are estimated by scaling
+* the variance on the weighted mean (rather than the weighted median)
+* by Pi/2. Otherwise, the output variances are calculated using the
+* COVEC array.
* ARR( NENT ) = DOUBLE PRECISION (Given)
* The list of ordered data for which the weighted median is
* required
@@ -47,7 +52,7 @@
* true in this array, otherwise the array is set to false.
* COVAR( * ) = DOUBLE PRECISION (Given)
* The packed variance-covariance matrix of the order statistics
-* from a normal distribution of size NENT.
+* from a normal distribution of size NENT. Not used if CALCMV is .FALSE.
* RESULT = DOUBLE PRECISION (Returned)
* The weighted median
* RESVAR = DOUBLE PRECISION (Returned)
@@ -77,6 +82,8 @@
* Added used array.
* 27-MAY-1992 (PDRAPER):
* Added direct variance estimates.
+* 27-AUG-2003 (DSB):
+* Added argument CALCMV.
* {enter_further_changes_here}
* Bugs:
@@ -91,6 +98,7 @@
INCLUDE 'SAE_PAR' ! Standard SAE constants
* Arguments Given:
+ LOGICAL CALCMV
INTEGER NENT
DOUBLE PRECISION ORDDAT( NENT )
DOUBLE PRECISION WEIGHT( NENT )
@@ -105,8 +113,13 @@
* Status:
INTEGER STATUS ! Global status
+* Local Constants:
+ DOUBLE PRECISION PIBY2 ! 90 degs in radians
+ PARAMETER ( PIBY2 = 1.57079632679489661923132 )
+
* Local Variables:
DOUBLE PRECISION TOTWT ! The total value of weights
+ DOUBLE PRECISION TRGWT ! The target value of weights
INTEGER I ! Loop variable
INTEGER J ! Loop variable
INTEGER K ! Loop variable
@@ -150,7 +163,7 @@
1 CONTINUE
* Search for median weight.
- TOTWT = TOTWT * 0.5D0
+ TRGWT = TOTWT * 0.5D0
WTSUM = 0.0D0
DO 2 I = 1, NENT
IF ( I .EQ. 1 ) THEN
@@ -160,7 +173,7 @@
: * 0.5D0
END IF
WTSUM = WTSUM + WTINC
- IF ( WTSUM .GT. TOTWT ) GO TO 66
+ IF ( WTSUM .GT. TRGWT ) GO TO 66
2 CONTINUE
66 I = MIN( I, NENT )
@@ -177,43 +190,51 @@
USED( I - 1 ) = .TRUE.
* Set weights factors
- W1 = ( WTSUM - TOTWT ) / MAX( WTINC, 1.0D-20 )
+ W1 = ( WTSUM - TRGWT ) / MAX( WTINC, 1.0D-20 )
W2 = 1.0D0 - W1
* Interpolate between data values
RESULT = D1 * W1 + D2 * W2
-* Sum the relevant ordered statistic variances and covariances.
-* weighting accordingly.
- VSUM = 0.0D0
- DO 3 K = LBND, UBND
- IF( K .EQ. LBND ) THEN
- IW = W1
- ELSE
- IW = W2
- END IF
- DO 4 J = K, UBND
- IF( J .EQ. LBND ) THEN
- JW = W1
+* If required, sum the relevant ordered statistic variances and
+* covariances weighting accordingly.
+ IF( CALCMV ) THEN
+ VSUM = 0.0D0
+ DO 3 K = LBND, UBND
+ IF( K .EQ. LBND ) THEN
+ IW = W1
ELSE
- JW = W2
+ IW = W2
END IF
+ DO 4 J = K, UBND
+ IF( J .EQ. LBND ) THEN
+ JW = W1
+ ELSE
+ JW = W2
+ END IF
* Sum variances and twice covariances ( off diagonal elements ).
- IF( K .EQ. J ) THEN
- VSUM = VSUM + IW * JW * COVAR( K + J * ( J - 1 )/ 2 )
- ELSE
- VSUM = VSUM +
- : 2.0D0 * IW * JW * COVAR( K + J * ( J - 1 )/ 2 )
- END IF
- 4 CONTINUE
- 3 CONTINUE
+ IF( K .EQ. J ) THEN
+ VSUM = VSUM + IW * JW * COVAR( K + J*( J - 1 )/2 )
+ ELSE
+ VSUM = VSUM +
+ : 2.0D0 * IW * JW * COVAR( K + J*( J - 1 )/2 )
+ END IF
+ 4 CONTINUE
+ 3 CONTINUE
* Right make the new variance estimate. Use the sum of variances
* and covariances of the order statistic of the `trimmed' sample size
* Sample variance changes to NENT * VAR to represent total variance
* of original data.
- RESVAR = VAR * NENT * VSUM
+ RESVAR = VAR * NENT * VSUM
+
+* If we are estimating the variance on the basis of the variance of the
+* mean, calculate the estimate.
+ ELSE
+ RESVAR = PIBY2/TOTWT
+ END IF
+
END IF
END
-* $Id: ccg1_wtm3.gen,v 1.1 1997/06/27 09:01:41 pwd Exp $
+* $Id: ccg1_wtm3d.f,v 1.3 2003/08/27 18:20:06 dsb Exp $
|
