supolar.c

1.Polarization analysis and filtering for three-component data 2.SUEIPOFI - EIgenimage (SVD) based

            /* write polarization attributes to files */
            if (rl) fputdata(rlfp, headerfp, data_rl, nt);
            if (theta) fputdata(thetafp, headerfp, data_theta, nt);
            if (phi) fputdata(phifp, headerfp, data_phi, nt);
            if (tau) fputdata(taufp, headerfp, data_tau, nt);
            if (ellip) {
                fputdata(e21fp, headerfp, data_e21, nt);
                fputdata(e31fp, headerfp, data_e31, nt);
                fputdata(e32fp, headerfp, data_e32, nt);
            }
            if (pln) fputdata(plnfp, headerfp, data_pln, nt);
            if (f1) fputdata(f1fp, headerfp, data_f1, nt);
            if (l1) fputdata(l1fp, headerfp, data_l1, nt);
            if (amp) {
                fputdata(erfp, headerfp, data_er, nt);
                fputdata(irfp, headerfp, data_ir, nt);
                fputdata(qrfp, headerfp, data_qr, nt);
            }
            if (dir) fputdata3c(dirfp, headerfp, data3c_dir, nt);
            
        } /* end of processing three-component dataset */
        
    } while (gettr(&tr)); /* end loop over traces */
    
    if (verbose) {
        fprintf(stderr,"\n");
        if (icomp) warn("last %d trace(s) skipped", icomp);
    }
    
    /* close files */
    efclose(headerfp);
    if (rl) efclose(rlfp);
    if (theta) efclose(thetafp);
    if (phi) efclose(phifp);
    if (tau) efclose(taufp);
    if (ellip) {
        efclose(e21fp);
        efclose(e31fp);
        efclose(e32fp);
    }
    if (pln) efclose(plnfp);
    if (f1) efclose(f1fp);
    if (l1) efclose(l1fp);
    if (amp) {
        efclose(erfp);
        efclose(irfp);
        efclose(qrfp);
    }
    if (dir) efclose(dirfp);
      
    return(CWP_Exit());
}



/**********************************************************************/
/* Functions used internally                                          */
/**********************************************************************/

/* calculate time window weights for a smooth time window,        */
/* after Sheriff, 1991, p. 335; tapered windows are suggested by  */
/* Jurkevics, 1988                                                */

void calc_window(float *w, int iwl, int iwin)
{
    int i,j;      /* loop indices within window */
    float m;      /* half of time window (iwl/2) */
    
    m = (float) iwl/2;
    
    switch (iwin) {
        case WBOXCAR :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 1.0;
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
        case WBARTLETT :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 1.0 - (fabs((float) i) / m);
                
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;
        case WHANNING :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = 0.5 + 0.5*cos(PI*fabs((float) i) / m);
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;            
        case WWELSH :
            for (i=-iwl/2,j=0;i<iwl/2;i++,j++) {
                w[j] = fabs( fabs((float) i) / m - 1.0);
                w[j] *= w[j];            
            }
            if (j<iwl) w[iwl-1] = w[0];
            break;            
    }
}


/* covariance of two time sequences *data1 and *data2,   */
/* calculated in a window of length iwl,                 */
/* weighted with factors *w and starting at index istart */

float covar(float *data1, float *data2, int istart, int iwl, float *w)
{
    register int i, j;
    float cov=0.0;
    float mean1=0.0;
    float mean2=0.0;
    
    for (i=istart,j=0;i<(istart+iwl);i++,j++) {
        mean1 += data1[i] * w[j];
        mean2 += data2[i] * w[j];
    }
    mean1 = mean1 / (float) iwl;
    mean2 = mean2 / (float) iwl;

    for (i=istart,j=0;i<(istart+iwl);i++,j++) {
        cov += (data1[i]-mean1) * (data2[i]-mean2) * w[j];
    }
    cov = cov / (float) iwl;
    return cov;
}   


/**********************************************************************/
/* Functions defining polarization attributes and angles              */
/**********************************************************************/

/* rectilinearity factor rl (3 different definitions) */

float calc_rl(float *d, float rlq, int opt)
{   
    float rl;
    
    if (!d[1]==0.0) {
        switch (opt) {
            default: 
                /* case 1: */
                /* rl definition after Kanasewich, 1981 */
                rl = 1 - pow(fabs(d[2]/d[1]), rlq);
                break;
            case 2:
                /* rl definition after Jurkevics, 1988 */
                rl = 1 - pow( 0.5*(fabs(d[2]/d[1]) + fabs(d[3]/d[1]) ), rlq);
                break;
            case 3:
                /* rl definition after Meyer, 1988 */
                rl = 1 - pow(fabs(d[2]/d[1]) + fabs(d[3]/d[1]), rlq);
                break;
        }
        return rl;
    }
    else {
        return 0.0;
    }
}

/* vertical (incidence) angle theta (different definitions) */
/* assumes a right handed coordinate system starting        */
/* with the vertical component Z (Z=1, N,R=2, E,T=3)        */

float calc_theta(float **v, int opt)
{
    float theta, horiz;
    
    switch (opt) {
        default:
            /* case 1: */
            /* definition after Kanasewich, 1981 */
            /* interval -pi/2 <= theta <= pi/2 */ 
            if (v[1][1]) {
                horiz = sqrt( v[2][1]*v[2][1] + v[3][1]*v[3][1] );
                theta = atan( horiz / v[1][1] );
            }
            else {
                theta = 0.0;
            }
            break;
        case 2:
        case 3:
            /* definition after Jurkevics, 1988 */
            /* interval 0.0 <= theta <= pi/2 */
            theta = acos( fabs(v[1][1]) );
            break;
    }
    return theta;            
}

/* horizontal azimuth phi (different definitions)    */
/* assumes a right handed coordinate system starting */
/* with the vertical component Z (Z=1, N,R=2, E,T=3) */

/* Note: The SIGN function is introduced to resolve  */
/* the 180 deg ambiguity by taking the positive      */
/* vertical component of v[1][1] (Jurkevics, 1988).  */
#define VSIGN ( (v[1][1]<0.0) ? -1.0 : 1.0 )

float calc_phi(float **v, int opt)
{
    float phi;
    
    switch (opt) {
        default:
            /* case 1: */
            /* definition after Kanasewich, 1981 */
            /* interval -pi/2 <= phi <= pi/2 */ 
            if (v[1][2]) {
                phi = atan( v[3][1] / v[2][1] );
            }
            else {
                phi = (v[3][1]>0.0) ? 0.5*PI : -0.5*PI;
            }
            break;
        case 2:
        case 3:
            /* definitions after Jurkevics, 1988 */
            /* interval -pi <= phi <= pi */
            if (v[2][1]) {
                phi = atan2( v[3][1]*VSIGN, v[2][1]*VSIGN);
            }
            else {
                phi = (v[3][1]>0.0) ? 0.5*PI*VSIGN : -0.5*PI*VSIGN;
            }
            
            /* interval 0.0 <= phi <= 2*pi */
            if (phi<0.0 && opt==3) phi += 2.0*PI;
            break;
    }   
    return phi;
}
#undef VSIGN

/* global polarization parameter tau (Samson, 1973) */

float calc_tau(float *d)
{
    float x1, x2, x3, x4, tau;
    
    if (d[1]) {   
        x1 = pow(( 1 - fabs(d[2]/d[1])), 2.);
        x2 = pow(( 1 - fabs(d[3]/d[1])), 2.);
        x3 = pow(( fabs(d[2]/d[1]) - fabs(d[3]/d[1]) ), 2.);
        x4 = pow(( 1 + fabs(d[2]/d[1]) - fabs(d[3]/d[1])), 2.);
        tau = sqrt( (x1+x2+x3) / (2*x4) );
        return tau;
    }
    else
        return 0.0;
}

/* ellipticities e_ik */

float calc_ellip(float *d, int d1, int d2)
{
    float ellip;
    
    if (d[1]) {
        ellip = sqrt( fabs(d[d1]/d[d2]) );
        return ellip;
    }
    else
        return 0.0;
}

/* planarity after Jurkevics, 1988 */

float calc_plan(float *d)
{
    float pln;
    
    if (d[1]+d[2]) {
        pln = 1.0 - 2.0*d[3] / (d[1] + d[2]);
        return pln;
    }
    else
        return 0.0;
}

/* flatness coefficient f1 after Benhama et. al, 1988 */

float calc_f1(float *d)
{
    float f1,x1,x2;
    
    x1 = 3.0 * calc_ellip(d,3,1);
    x2 = 1.0 + calc_ellip(d,2,1) + calc_ellip(d,3,1);
    f1 = 1.0 - x1 / x2; 
    return f1;
}

/* linearity coefficient l1 */

float calc_l1(float *d)
{
    float l1,x1,x2;
    
    x1 = 3. * ( calc_ellip(d,2,1) + calc_ellip(d,3,1) );
    x2 = 2. * ( 1 + calc_ellip(d,2,1) + calc_ellip(d,3,1));
    l1 = 1. - x1 / x2;
    return l1;
}

/* direction cosines or directivity functions (3 components) */

void calc_dir(float **data3c_dir, float **v, int it)
{
    int i;
    for (i=1;i<=3;i++) {
        data3c_dir[i][it]=v[i][1];
    }
}

/* amplitude parameters */

/* eigenresultant */

float calc_er(float *d)
{
    float er;
    er = sqrt(fabs(d[1]));
    return er;
}

/* instantaneous and quadratic resultant (Meyer, 1988) */

void ampparams(float **indata, float *data_ir, float *data_qr, int nt, int iwl)
{
    int i,it;
    float sqrsum;
    
    for (it=0;it<nt;it++) {
        /* instantaneous resultant */
        data_ir[it] = sqrt(indata[1][it]*indata[1][it]+ \
            indata[2][it]*indata[2][it]+indata[3][it]*indata[3][it]);
        
        /* quadratic resultant */
        if ( (it >= iwl/2) && (it < nt-iwl/2) ) {
            sqrsum=0.0;
            for (i=it-iwl/2;i<(it+iwl);i++) {
                sqrsum += indata[1][i]*indata[1][i]+ \
                    indata[2][i]*indata[2][i]+indata[3][i]*indata[3][i];
            }
            data_qr[it] = sqrsum / (float) iwl;
        }
    }
}


/**********************************************************************/
/* Functions for data output                                          */
/**********************************************************************/


/* write one-component data into file */

void fputdata(FILE *fileptr, FILE *headerptr, float *outdata, int nt)
{    
    efread(&tr, 1, HDRBYTES, headerptr);
    erewind(headerptr);
    
    memcpy((void *)tr.data, (const void *) outdata, nt*FSIZE);

    fputtr(fileptr, &tr);
}

/* write three-component data into file */

void fputdata3c(FILE *fileptr, FILE *headerptr, float **outdata3c, int nt)
{
    int i;
   
    for(i=1;i<=3;i++) {
        efread(&tr, 1, HDRBYTES, headerptr);
        
        memcpy((void *)tr.data, (const void *) outdata3c[i], nt*FSIZE);

        fputtr(fileptr, &tr);
    }
    erewind(headerptr);
}

/* END OF FILE */