00004 {
00005    static char sval[128] ;
00006    double aval=fabs(val) , tval ;
00007    int    lv , qv ;
00008 
00009    if( aval == 0.0 ){ strcpy(sval,"Zero"); return; }
00010 
00011    if( val < 0.0 ){ strcpy(sval,"-"); } else { sval[0] = '\0'; }
00012 
00013    lv   = (int) floor(log10(aval)/3.0) ;
00014    tval = pow(10.0,(double)(3*lv)) ;
00015    qv   = (int) rint(aval/tval) ;
00016    sprintf( sval+strlen(sval) , "%d" , qv ) ;
00017 
00018    switch( lv ){
00019 
00020      case 0: break ;
00021 
00022      case 1: strcat(sval+strlen(sval)," thousand")    ; break ;
00023      case 2: strcat(sval+strlen(sval)," million" )    ; break ;
00024      case 3: strcat(sval+strlen(sval)," billion" )    ; break ;
00025      case 4: strcat(sval+strlen(sval)," trillion")    ; break ;
00026      case 5: strcat(sval+strlen(sval)," quadrillion") ; break ;
00027      case 6: strcat(sval+strlen(sval)," quintillion") ; break ;
00028 
00029      case -1: strcat(sval+strlen(sval)," thousand-ths") ; break ;
00030      case -2: strcat(sval+strlen(sval)," million-ths")  ; break ;
00031      case -3: strcat(sval+strlen(sval)," billion-ths")  ; break ;
00032      case -4: strcat(sval+strlen(sval)," trillion-ths") ; break ;
00033 
00034      default:
00035        strcat(sval+strlen(sval)," jillion") ;
00036        if( lv < 0 ) strcat(sval+strlen(sval),"-ths") ;
00037      break ;
00038    }
00039 
00040    return (char *)sval ;
00041 }

00038 {
00039    /* loop counters */
00040 
00041    int jj , ii ;
00042 
00043    /* variables for CL1 (types declared in f2c.h) */
00044 
00045    integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00046    real *q , toler , *x , *res , error , *cu ;
00047 
00048    /*-- check inputs --*/
00049 
00050    if( ndim < 1 || nvec < 1 )                         kode = 4 ;
00051    if( A == NULL || y == NULL || z == NULL )          kode = 4 ;
00052    for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00053 
00054    if( kode ){
00055      fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00056      return (float)(-kode) ;
00057    }
00058 
00059    /*-- setup call to CL1 --*/
00060 
00061    k     = ndim ;
00062    l     = 0 ;     /* no linear equality constraints */
00063    m     = 0 ;     /* no linear inequality constraints */
00064    n     = nvec ;
00065 
00066    klmd  = k+l+m ;
00067    klm2d = k+l+m+2 ;
00068    nklmd = n+k+l+m ;
00069    n2d   = n+2 ;
00070 
00071    kode  = (cony != 0) ; /* enforce implicit constraints on x[] */
00072    iter  = 11*klmd ;
00073 
00074    toler = 0.0001 ;
00075    error = 0.0 ;
00076 
00077    /* input/output matrices & vectors */
00078 
00079    q     = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00080    x     = (real *) calloc( 1, sizeof(real) * n2d ) ;
00081    res   = (real *) calloc( 1, sizeof(real) * klmd ) ;
00082 
00083    /* workspaces */
00084 
00085    cu    = (real *)    calloc( 1, sizeof(real) * 2*nklmd ) ;
00086    iu    = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00087    s     = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00088 
00089    /* load matrices & vectors */
00090 
00091    for( jj=0 ; jj < nvec ; jj++ )
00092       for( ii=0 ; ii < ndim ; ii++ )
00093          q[ii+jj*klm2d] = A[jj][ii] ;   /* matrix */
00094 
00095    for( ii=0 ; ii < ndim ; ii++ )
00096       q[ii+nvec*klm2d] = z[ii] ;        /* vector */
00097 
00098    if( cony ){
00099      for( jj=0 ; jj < nvec ; jj++ )     /* constraints on solution */
00100        x[jj] = (y[jj] < 0.0) ? -1.0
00101               :(y[jj] > 0.0) ?  1.0 : 0.0 ;
00102    }
00103 
00104    for( ii=0 ; ii < klmd ; ii++ )       /* no constraints on resids */
00105       res[ii] = 0.0 ;
00106 
00107    /*-- do the work --*/
00108 
00109    cl1_fort( &k, &l, &m, &n,
00110              &klmd, &klm2d, &nklmd, &n2d,
00111              q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00112 
00113    free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
00114 
00115    if( kode != 0 ){
00116      free(x) ;
00117 #if 0
00118      switch( kode ){
00119        case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00120        case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n")     ; break;
00121        case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n")      ; break;
00122       default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n")     ; break;
00123      }
00124 #endif
00125      return (float)(-kode) ;
00126    }
00127 
00128    /*-- copy results into output --*/
00129 
00130    for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00131 
00132    free(x) ; return (float)error ;
00133 }

00172 {
00173    /* loop counters */
00174 
00175    int jj , ii ;
00176 
00177    /* variables for CL1 (types declared in f2c.h) */
00178 
00179    integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00180    real *q , toler , *x , *res , error , *cu ;
00181 
00182    /*-- check inputs --*/
00183 
00184    if( ndim < 1 || nvec < 1 )                         kode = 4 ;
00185    if( A == NULL || y == NULL || z == NULL )          kode = 4 ;
00186    for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00187 
00188    if( kode ){
00189      fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00190      return (float)(-kode) ;
00191    }
00192 
00193    /*-- setup call to CL1 --*/
00194 
00195    k     = ndim ;
00196    l     = 0 ;     /* no linear equality constraints */
00197    m     = 0 ;     /* no linear inequality constraints */
00198    n     = nvec ;
00199 
00200    klmd  = k+l+m ;
00201    klm2d = k+l+m+2 ;
00202    nklmd = n+k+l+m ;
00203    n2d   = n+2 ;
00204 
00205    kode  = (cony != 0 || conr != 0) ; /* enforce implicit constraints on x[] */
00206    iter  = 11*klmd ;
00207 
00208    toler = 0.0001 ;
00209    error = 0.0 ;
00210 
00211    /* input/output matrices & vectors */
00212 
00213    q     = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00214    x     = (real *) calloc( 1, sizeof(real) * n2d ) ;
00215    res   = (real *) calloc( 1, sizeof(real) * klmd ) ;
00216 
00217    /* workspaces */
00218 
00219    cu    = (real *)    calloc( 1, sizeof(real) * 2*nklmd ) ;
00220    iu    = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00221    s     = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00222 
00223    /* load matrices & vectors */
00224 
00225    for( jj=0 ; jj < nvec ; jj++ )
00226       for( ii=0 ; ii < ndim ; ii++ )
00227          q[ii+jj*klm2d] = A[jj][ii] ;   /* matrix */
00228 
00229    for( ii=0 ; ii < ndim ; ii++ )
00230       q[ii+nvec*klm2d] = z[ii] ;        /* vector */
00231 
00232    if( cony ){
00233      for( jj=0 ; jj < nvec ; jj++ )     /* constraints on solution */
00234        x[jj] = (y[jj] < 0.0) ? -1.0
00235               :(y[jj] > 0.0) ?  1.0 : 0.0 ;
00236    }
00237 
00238    if( conr ){
00239      for( ii=0 ; ii < ndim ; ii++ )     /* constraints on resids */
00240        res[ii] = (rez[ii] < 0.0) ? -1.0
00241                 :(rez[ii] > 0.0) ?  1.0 : 0.0 ;
00242    }
00243 
00244    /*-- do the work --*/
00245 
00246    cl1_fort( &k, &l, &m, &n,
00247              &klmd, &klm2d, &nklmd, &n2d,
00248              q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00249 
00250    free(q) ; free(cu) ; free(iu) ; free(s) ;
00251 
00252    if( kode != 0 ){
00253      free(x) ; free(res) ;
00254 #if 0
00255      switch( kode ){
00256        case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00257        case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n")     ; break;
00258        case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n")      ; break;
00259       default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n")     ; break;
00260      }
00261 #endif
00262      return (float)(-kode) ;
00263    }
00264 
00265    /*-- copy results into output --*/
00266 
00267    for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00268 
00269    for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = (float) res[ii] ;
00270 
00271    free(res); free(x); return (float)error;
00272 }

00302 {
00303    if( xx == NULL || ww == NULL ) return ;
00304    if( n < 2 || n > 20  ){ *xx = *ww = NULL ; return ; }
00305 
00306    *xx = xlag[n] ; *ww = wlag[n] ; return ;
00307 }

00032 {
00033    int ii,jj , nfac , *fac ;
00034    int fd ; FILE *fp ;
00035    char qbuf[128] ;
00036 
00037 #ifndef DONT_USE_MKSTEMP
00038    char fname[] = "/tmp/afniXXXXXX" ;
00039 #else
00040    char *fname ;
00041 #endif
00042 
00043    if( npt < 3 || xyz == NULL || ijk == NULL ){
00044       fprintf(stderr,"qhull_wrap: bad inputs\n") ;
00045       return 0 ;
00046    }
00047 
00048 #ifndef DONT_USE_MKSTEMP
00049    fd = mkstemp( fname ) ;
00050    if( fd == -1 ){ fprintf(stderr,"qhull_wrap: mkstemp fails\n"); return 0; }
00051    fp = fdopen( fd , "w" ) ;
00052    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: fdopen fails\n"); close(fd); return 0; }
00053 #else
00054    fname = tmpnam(NULL) ;
00055    if( fname == NULL ){ fprintf(stderr,"qhull_wrap: tmpnam fails\n"); return 0; }
00056    fp = fopen( fname , "w" ) ;
00057    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: fopen fails\n"); return 0; }
00058 #endif
00059 
00060    fprintf(fp,"3\n%d\n",npt) ;
00061    for( ii=0 ; ii < npt ; ii++ )
00062       fprintf(fp,"%g %g %g\n",xyz[3*ii],xyz[3*ii+1],xyz[3*ii+2]) ;
00063 
00064    fclose(fp) ;
00065 
00066    sprintf(qbuf,"qhull -i -Pp < %s",fname) ;
00067    fp = popen( qbuf , "r" ) ;
00068    if( fp == NULL ){ fprintf(stderr,"qhull_wrap: popen fails\n"); remove(fname); return 0; }
00069 
00070    jj = fscanf(fp,"%d",&nfac) ;
00071    if( jj != 1 || nfac < 1 ){ fprintf(stderr,"qhull_wrap: 1st fscanf fails\n"); pclose(fp); remove(fname); return 0; }
00072 
00073    fac = (int *) malloc( sizeof(int)*3*nfac ) ;
00074    if( fac == NULL ){ fprintf(stderr,"qhull_wrap: malloc fails\n"); pclose(fp); remove(fname); return 0; }
00075 
00076    for( ii=0 ; ii < nfac ; ii++ ){
00077       jj = fscanf(fp,"%d %d %d",fac+(3*ii),fac+(3*ii+1),fac+(3*ii+2)) ;
00078       if( jj < 3 ){
00079          fprintf(stderr,"qhull_wrap: fscanf fails at ii=%d\n",ii) ;
00080          pclose(fp); remove(fname); free(fac); return 0;
00081       }
00082    }
00083 
00084    pclose(fp); remove(fname);
00085 
00086    *ijk = fac ; return nfac ;
00087 }

00033 {
00034    register int i , j ;           /* scanning indices */
00035    register float temp , pivot ;  /* holding places */
00036    register float * a = ar ;
00037 
00038    int left , right , mid , nodd ;
00039 
00040    switch( n ){
00041       case 1: return ar[0] ;
00042       case 2: return 0.5*(ar[0]+ar[1]) ;
00043       case 3: return MED3( ar[0] , ar[1] , ar[2] ) ;
00044       case 4: return median_float4( ar[0],ar[1],ar[2],ar[3] ) ;
00045       case 5: return median_float5( ar ) ;
00046       case 7: return median_float7( ar ) ;
00047       case 9: return median_float9( ar ) ;
00048    }
00049 
00050    left = 0 ; right = n-1 ; mid = n/2 ; nodd = ((n & 1) != 0) ;
00051 
00052    /* loop while the subarray is at least 3 long */
00053 
00054    while( right-left > 1  ){  /* work on subarray from left -> right */
00055 
00056       i = ( left + right ) / 2 ;   /* middle of subarray */
00057 
00058       /* sort the left, middle, and right a[]'s */
00059 
00060       if( a[left] > a[i]     ) SWAP( a[left]  , a[i]     ) ;
00061       if( a[left] > a[right] ) SWAP( a[left]  , a[right] ) ;
00062       if( a[i] > a[right]    ) SWAP( a[right] , a[i]     ) ;
00063 
00064       pivot = a[i] ;                 /* a[i] is the median-of-3 pivot! */
00065       a[i]  = a[right] ;
00066 
00067       i = left ;                     /* initialize scanning */
00068       j = right ;
00069 
00070       /*----- partition:  move elements bigger than pivot up and elements
00071                           smaller than pivot down, scanning in from ends -----*/
00072 
00073       do{
00074         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00075         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00076 
00077         if( j <= i ) break ;         /* if j meets i, quit */
00078 
00079         SWAP( a[i] , a[j] ) ;
00080       } while( 1 ) ;
00081 
00082       /*----- at this point, the array is partitioned -----*/
00083 
00084       a[right] = a[i] ;           /* restore the pivot */
00085       a[i]     = pivot ;
00086 
00087       if( i == mid ){             /* good luck */
00088          if( nodd ) return pivot ;   /* exact middle of array */
00089 
00090          temp = a[left] ;            /* must find next largest element below */
00091          for( j=left+1 ; j < i ; j++ )
00092             if( a[j] > temp ) temp = a[j] ;
00093 
00094          return 0.5*(pivot+temp) ;   /* return average */
00095       }
00096 
00097       if( i <  mid ) left  = i ; /* throw away bottom partition */
00098       else           right = i ; /* throw away top partition    */
00099 
00100    }  /* end of while sub-array is long */
00101 
00102    return (nodd) ? a[mid] : 0.5*(a[mid]+a[mid-1]) ;
00103 }

00110 {
00111    float * q = (float *) malloc(sizeof(float)*n) ;
00112    float me,ma ;
00113    register int ii ;
00114 
00115    memcpy(q,ar,sizeof(float)*n) ;  /* duplicate input array */
00116    me = qmed_float( n , q ) ;      /* compute median */
00117 
00118    for( ii=0 ; ii < n ; ii++ )     /* subtract off median */
00119       q[ii] = fabs(q[ii]-me) ;     /* (absolute deviation) */
00120 
00121    ma = qmed_float( n , q ) ;      /* MAD = median absolute deviation */
00122 
00123    free(q) ;                       /* 05 Nov 2001 */
00124 
00125    *med = me ; *mad = ma ; return ;
00126 }

00134 {
00135    qsrec_doubleint( n , a , ia , QS_CUTOFF ) ;
00136    isort_doubleint( n , a , ia ) ;
00137    return ;
00138 }

00135 {
00136    qsrec_floatfloat( n , a , ia , QS_CUTOFF ) ;
00137    isort_floatfloat( n , a , ia ) ;
00138    return ;
00139 }

00134 {
00135    qsrec_floatint( n , a , ia , QS_CUTOFF ) ;
00136    isort_floatint( n , a , ia ) ;
00137    return ;
00138 }

00137 {
00138    qsrec_floatstuff( n , a , ia , QS_CUTOFF ) ;
00139    isort_floatstuff( n , a , ia ) ;
00140    return ;
00141 }

00134 {
00135    qsrec_intint( n , a , ia , QS_CUTOFF ) ;
00136    isort_intint( n , a , ia ) ;
00137    return ;
00138 }

00102 {
00103    float *xyz ;
00104    double cp,sp,ca,sa ;
00105    int ii ;
00106 
00107    if( npt < 3 || pol == NULL || azi == NULL || wt == NULL ){
00108       fprintf(stderr,"sphere_voronoi_angles: bad inputs\n"); return 0;
00109    }
00110 
00111    /* make 3-vectors of the points on the sphere */
00112 
00113    xyz = (float *) malloc( sizeof(float) * (3*npt) ) ;
00114 
00115    for( ii=0 ; ii < npt ; ii++ ){
00116       cp = cos(pol[ii]) ; sp = sin(pol[ii]) ;
00117       ca = cos(azi[ii]) ; sa = sin(azi[ii]) ;
00118       xyz[3*ii]   = sp * ca ;
00119       xyz[3*ii+1] = sp * sa ;
00120       xyz[3*ii+2] = cp ;
00121    }
00122 
00123    ii = sphere_voronoi_vectors( npt , xyz , wt ) ;
00124    free(xyz) ; return ii ;
00125 }

00135 {
00136    int ntri , *tri , ii,jj , pp,qq,rr ;
00137    float *ww ;
00138    double cp,sp,ca,sa ;
00139    double xpq,ypq,zpq , xpr,ypr,zpr , xqr,yqr,zqr , xcc,ycc,zcc ;
00140    double xpp,ypp,zpp , xqq,yqq,zqq , xrr,yrr,zrr , xnn,ynn,znn ;
00141    double pp_pq , pp_pr , pp_cc ,
00142           qq_pq , qq_qr , qq_cc ,
00143           rr_qr , rr_cc , rr_pr ,
00144           pq_cc , qr_cc , pr_cc , ss ;
00145    double a_pp_pq_cc , a_pp_pr_cc ,
00146           a_qq_pq_cc , a_qq_qr_cc ,
00147           a_rr_qr_cc , a_rr_pr_cc  ;
00148 
00149    if( npt < 3 || xyz == NULL || wt == NULL ){
00150       fprintf(stderr,"sphere_voronoi: bad inputs\n"); return 0;
00151    }
00152 
00153    /* compute convex hull triangular facets */
00154 
00155    ntri = qhull_wrap( npt , xyz , &tri ) ;
00156    if( ntri == 0 ){ fprintf(stderr,"sphere_voronoi: qhull_wrap fails\n"); free(xyz); return 0; }
00157 
00158    /* initialize output */
00159 
00160    ww = (float *) malloc( sizeof(float) * npt ) ;
00161    for( ii=0 ; ii < npt ; ii++ ) ww[ii] = 0.0 ;
00162 
00163    for( jj=0 ; jj < ntri ; jj++ ){  /* loop over triangles */
00164 
00165 
00166       /* triangle vertices pp,qq,rr      * pp
00167                                         / \
00168                                        /   \
00169                                    pq *     * pr
00170                                      /   *   \
00171                                     /   cc    \
00172                                 qq *-----------* rr   */
00173 
00174       pp  = tri[3*jj  ] ;
00175       xpp = xyz[3*pp  ] ; ypp = xyz[3*pp+1] ; zpp = xyz[3*pp+2] ;
00176       qq  = tri[3*jj+1] ;
00177       xqq = xyz[3*qq  ] ; yqq = xyz[3*qq+1] ; zqq = xyz[3*qq+2] ;
00178       rr  = tri[3*jj+2] ;
00179       xrr = xyz[3*rr  ] ; yrr = xyz[3*rr+1] ; zrr = xyz[3*rr+2] ;
00180 
00181       /* midpoints pq,pr,qr, and centroid cc */
00182       /*** Q: should centroid be replaced by normal?
00183               what about orientation, largeness?     ***/
00184 
00185       xpq = 0.5*(xpp+xqq) ; ypq = 0.5*(ypp+yqq) ; zpq = 0.5*(zpp+zqq) ;
00186       xpr = 0.5*(xpp+xrr) ; ypr = 0.5*(ypp+yrr) ; zpr = 0.5*(zpp+zrr) ;
00187       xqr = 0.5*(xqq+xrr) ; yqr = 0.5*(yqq+yrr) ; zqr = 0.5*(zqq+zrr) ;
00188 
00189       xcc = 0.3333333*(xpp+xqq+xrr) ;
00190       ycc = 0.3333333*(ypp+yqq+yrr) ;
00191       zcc = 0.3333333*(zpp+zqq+zrr) ;
00192 
00193 #undef SCL
00194 #define SCL(a,b,c) 1.0/sqrt(a*a+b*b+c*c)
00195 
00196 #undef USE_NORMAL
00197 #ifdef USE_NORMAL
00198 # define XCROSS(a,b,c,x,y,z) ((b)*(z)-(c)*(y))
00199 # define YCROSS(a,b,c,x,y,z) ((c)*(x)-(a)*(z))
00200 # define ZCROSS(a,b,c,x,y,z) ((a)*(y)-(b)*(x))
00201       { double apq=xpp-xqq , bpq=ypp-yqq , cpq=zpp-zqq ,
00202                aqr=xqq-xrr , bqr=yqq-yrr , cqr=zqq-zrr  ;
00203 
00204         xnn = XCROSS(apq,bpq,cpq,aqr,bqr,cqr) ,
00205         ynn = YCROSS(apq,bpq,cpq,aqr,bqr,cqr) ,
00206         znn = ZCROSS(apq,bpq,cpq,aqr,bqr,cqr)  ;
00207 
00208         cp = SCL(xnn,ynn,znn) ; xnn *= cp ; ynn *= cp ; znn *= cp ;
00209         if( xnn*xcc + ynn*ycc + znn*zcc < 0 ){
00210            xnn = -xnn ; ynn = -ynn ; znn = -znn ;
00211         }
00212       }
00213 
00214 # define xVV xnn
00215 # define yVV ynn
00216 # define zVV znn
00217 
00218 #else
00219 
00220 # define xVV xcc
00221 # define yVV ycc
00222 # define zVV zcc
00223 
00224 #endif
00225 
00226       /* project pq,pr,qr,cc to sphere (nn is already on sphere) */
00227 
00228       cp = SCL(xpq,ypq,zpq) ; xpq *= cp ; ypq *= cp ; zpq *= cp ;
00229       cp = SCL(xpr,ypr,zpr) ; xpr *= cp ; ypr *= cp ; zpr *= cp ;
00230       cp = SCL(xqr,yqr,zqr) ; xqr *= cp ; yqr *= cp ; zqr *= cp ;
00231       cp = SCL(xcc,ycc,zcc) ; xcc *= cp ; ycc *= cp ; zcc *= cp ;
00232 
00233 #undef  ANG
00234 #define ANG(u1,u2,u3,v1,v2,v3) acos(u1*v1+u2*v2+u3*v3)
00235 
00236       /* compute distance on surface between points:
00237          aa_bb = between points aa and bb, from the picture above */
00238 
00239       pp_pq = ANG( xpp,ypp,zpp , xpq,ypq,zpq ) ;
00240       pp_cc = ANG( xpp,ypp,zpp , xVV,yVV,zVV ) ;
00241       pp_pr = ANG( xpp,ypp,zpp , xpr,ypr,zpr ) ;
00242 
00243       qq_pq = ANG( xqq,yqq,zqq , xpq,ypq,zpq ) ;
00244       qq_qr = ANG( xqq,yqq,zqq , xqr,yqr,zqr ) ;
00245       qq_cc = ANG( xqq,yqq,zqq , xVV,yVV,zVV ) ;
00246 
00247       rr_qr = ANG( xrr,yrr,zrr , xqr,yqr,zqr ) ;
00248       rr_pr = ANG( xrr,yrr,zrr , xpr,ypr,zpr ) ;
00249       rr_cc = ANG( xrr,yrr,zrr , xVV,yVV,zVV ) ;
00250 
00251       pq_cc = ANG( xpq,ypq,zpq , xVV,yVV,zVV ) ;
00252       qr_cc = ANG( xqr,yqr,zqr , xVV,yVV,zVV ) ;
00253       pr_cc = ANG( xpr,ypr,zpr , xVV,yVV,zVV ) ;
00254 
00255       /* for each vertex,
00256          compute areas of 2 subtriangles it touches,
00257          add these into the area weight for that vertex */
00258 
00259 #undef  SS
00260 #define SS(a,b,c) ((a+b+c)/2)
00261 #undef  ATR
00262 #define ATR(s,a,b,c) (4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2))))
00263 
00264       ss      = SS(pp_pq,pp_cc,pq_cc) ;     /* subtriangle pp,pq,cc */
00265       ww[pp] += a_pp_pq_cc = ATR(ss,pp_pq,pp_cc,pq_cc) ;
00266 
00267       ss      = SS(pp_pr,pp_cc,pr_cc) ;     /* subtriangle pp,pr,cc */
00268       ww[pp] += a_pp_pr_cc = ATR(ss,pp_pr,pp_cc,pr_cc) ;
00269 
00270       ss      = SS(qq_pq,qq_cc,pq_cc) ;     /* subtriangle qq,pq,cc */
00271       ww[qq] += a_qq_pq_cc = ATR(ss,qq_pq,qq_cc,pq_cc) ;
00272 
00273       ss      = SS(qq_qr,qq_cc,qr_cc) ;     /* subtriangle qq,qr,cc */
00274       ww[qq] += a_qq_qr_cc = ATR(ss,qq_qr,qq_cc,qr_cc) ;
00275 
00276       ss      = SS(rr_qr,rr_cc,qr_cc) ;     /* subtriangle rr,qr,cc */
00277       ww[rr] += a_rr_qr_cc = ATR(ss,rr_qr,rr_cc,qr_cc) ;
00278 
00279       ss      = SS(rr_pr,rr_cc,pr_cc) ;     /* subtriangle rr,pr,cc */
00280       ww[rr] += a_rr_pr_cc = ATR(ss,rr_pr,rr_cc,pr_cc) ;
00281 
00282 
00283 #if 0          /* debugging printouts */
00284 # undef  DDD
00285 # define DDD(x,y,z,a,b,c) sqrt((x-a)*(x-a)+(y-b)*(y-b)+(z-c)*(z-c))
00286 
00287       cp = DDD(xpp,ypp,zpp,xqq,yqq,zqq) ;
00288       sp = DDD(xpp,ypp,zpp,xrr,yrr,zrr) ;
00289       ca = DDD(xqq,yqq,zqq,xrr,yrr,zrr) ;
00290       sa = (cp+sp+ca)/2 ;
00291       ss = sqrt(sa*(sa-cp)*(sa-sp)*(sa-ca)) ;
00292 
00293       fprintf(stderr,"triangle %d: pp=%d qq=%d rr=%d  AREA=%6.3f PLANAR=%6.3f\n"
00294                      "  xpp=%6.3f ypp=%6.3f zpp=%6.3f\n"
00295                      "  xqq=%6.3f yqq=%6.3f zqq=%6.3f\n"
00296                      "  xrr=%6.3f yrr=%6.3f zrr=%6.3f\n"
00297                      "  xpq=%6.3f ypq=%6.3f zpq=%6.3f\n"
00298                      "  xqr=%6.3f yqr=%6.3f zqr=%6.3f\n"
00299                      "  xpr=%6.3f ypr=%6.3f zpr=%6.3f\n"
00300                      "  xcc=%6.3f ycc=%6.3f zcc=%6.3f\n"
00301 #ifdef USE_NORMAL
00302                      "  xnn=%6.3f ynn=%6.3f znn=%6.3f\n"
00303 #endif
00304                      "  pp_pq=%6.3f pp_pr=%6.3f pp_cc=%6.3f\n"
00305                      "  qq_pq=%6.3f qq_qr=%6.3f qq_cc=%6.3f\n"
00306                      "  rr_qr=%6.3f rr_cc=%6.3f rr_pr=%6.3f\n"
00307                      "  pq_cc=%6.3f qr_cc=%6.3f pr_cc=%6.3f\n"
00308                      "  a_pp_pq_cc=%6.3f a_pp_pr_cc=%6.3f\n"
00309                      "  a_qq_pq_cc=%6.3f a_qq_qr_cc=%6.3f\n"
00310                      "  a_rr_qr_cc=%6.3f a_rr_pr_cc=%6.3f\n" ,
00311                jj,pp,qq,rr ,
00312                a_pp_pq_cc+a_pp_pr_cc+a_qq_pq_cc+a_qq_qr_cc+a_rr_qr_cc+a_rr_pr_cc , ss ,
00313                xpp, ypp, zpp,
00314                xqq, yqq, zqq,
00315                xrr, yrr, zrr,
00316                xpq, ypq, zpq,
00317                xqr, yqr, zqr,
00318                xpr, ypr, zpr,
00319                xcc, ycc, zcc,
00320 #ifdef USE_NORMAL
00321                xnn, ynn, znn,
00322 #endif
00323                pp_pq, pp_pr, pp_cc,
00324                qq_pq, qq_qr, qq_cc,
00325                rr_qr, rr_cc, rr_pr,
00326                pq_cc, qr_cc, pr_cc,
00327                a_pp_pq_cc, a_pp_pr_cc,
00328                a_qq_pq_cc, a_qq_qr_cc,
00329                a_rr_qr_cc, a_rr_pr_cc  ) ;
00330 #endif
00331 
00332    } /* end of loop over triangles */
00333 
00334    /* exit stage left */
00335 
00336    *wt = ww ; return 1 ;
00337 }

00082 {
00083    integer mm,nn , lda,ldu,ldv , ierr ;
00084    doublereal *aa, *ww , *uu , *vv , *rv1 ;
00085    logical    matu , matv ;
00086 
00087    if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
00088 
00089    mm  = m ;
00090    nn  = n ;
00091    aa  = a ;
00092    lda = m ;
00093    ww  = s ;
00094 
00095    if( u == NULL ){
00096      matu = (logical) 0 ;
00097      uu   = (doublereal *)malloc(sizeof(double)*m*n) ;
00098    } else {
00099      matu = (logical) 1 ;
00100      uu = u ;
00101    }
00102    ldu = m ;
00103 
00104    if( v == NULL ){
00105      matv = (logical) 0 ;
00106      vv   = NULL ;
00107    } else {
00108      matv = (logical) 1 ;
00109      vv = v ;
00110    }
00111    ldv = n ;
00112 
00113    rv1 = (double *) malloc(sizeof(double)*n) ;  /* workspace */
00114 
00115    (void) svd_( &mm , &nn , &lda , aa , ww ,
00116                 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
00117 
00118    free((void *)rv1) ;
00119 
00120    if( u == NULL ) free((void *)uu) ;
00121    return ;
00122 }

00026 {
00027    integer nm , matz , ierr ;
00028    double *fv1 , *fv2 ;
00029 
00030    if( a == NULL || e == NULL || n < 1 ) return ;
00031 
00032    if( n == 1 ){
00033      e[0] = a[0] ; a[0] = 1.0 ; return ;  /* degenerate case */
00034    }
00035 
00036    fv1 = (double *) malloc(sizeof(double)*n) ;  /* workspaces */
00037    fv2 = (double *) malloc(sizeof(double)*n) ;
00038 
00039    nm = n ; matz = 1 ; ierr = 0 ;
00040 
00041    rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
00042 
00043    free((void *)fv1) ; free((void *)fv2) ;
00044    return ;
00045 }

00050 {
00051    integer nm , matz , ierr ;
00052    double *fv1 , *fv2 ;
00053 
00054    if( a == NULL || e == NULL || n < 1 ) return ;
00055 
00056    if( n == 1 ){ e[0] = a[0] ; return ; } /* degenerate case */
00057 
00058    fv1 = (double *) malloc(sizeof(double)*n) ;  /* workspaces */
00059    fv2 = (double *) malloc(sizeof(double)*n) ;
00060 
00061    nm = n ; matz = 0 ; ierr = 0 ;
00062 
00063    rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
00064 
00065    free((void *)fv1) ; free((void *)fv2) ;
00066    return ;
00067 }