Thread: An update

Started: 2022-07-29 21:46:12
Last activity: 2022-07-29 21:46:12
Topics: SAC Help
刘有山
2022-07-29 21:46:12
Dear Mrs./Miss,

Thank you for your reading.

I have fixed those bugs mentioned in the last Email.

In addition, I have added a parameter for xtransfer.c and transfer.c. As a result, a water_level parameter is allowed from user.

Using the water_level, I compute threshold = h_max * pow(10., -water_level/20.0), which is similar to those done in obspy.

The default water_level is set as 600.

The involved files are as follows, /inc/icm.h, /src/icm/transfer.c, /src/icm/xtransfer.c.

I have uploaded the updated files in the attachments.

If you accept it, you can replace them.

Thanks a lot!

Best Wishes,

Youshan Liu
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <float.h>

#include "amf.h"
#include "icm.h"
#include "hdr.h"
#include "co.h"
#include "EVRESPnames.h"

#include "cssListOps/dblPublicDefs.h"
#include "smDataIO.h"
#include "cssListOps/cssStrucs.h"
#define EPS 1e-6

#define EQUAL(x,y) ( fabs(x - y) < EPS )

int GetCalibCalperFromWfdisc(DBlist tree, int wfid, double *calib,
double *calper);

int
NearestInt(double D) {
double Floor = floor(D);
double Ceil = ceil(D);
return (Ceil - D > D - Floor) ? (int) Floor : (int) Ceil;
}

int
IsNormalized(double calper, int nfreq, double delfrq, const double *xre,
const double *xim) {
double value;
double NormFreq = 1.0 / calper;
int j = NearestInt(NormFreq / delfrq);
if (j < 0 || j >= nfreq) {
printf
("\tCannot determine normalization status because NormFreq is outside bounds of transfer function.\n");
return 0;
}
value = sqrt(xre[j] * xre[j] + xim[j] * xim[j]);
return fabs(value - 1.0) < 0.01;
}

static int
is_type(char *name, char *type) {
return(strncasecmp(name, type, strlen(type)) == 0);
}

double
GetNormalizationFactor(int nfreq, double delfrq, const double *xre,
const double *xim) {
/* This function attempts to determine a scaling factor to apply
to the seismogram data to account for normalized vs non-normalized
transfer functions, and scaled vs unscaled data.
When the transfer function is non-normalized and the data were scaled
on input, the scale factor is 1/calib.
When the transfer function is normalized and the data have not been
scaled on input, the scale factor is calib.
In all other cases, the scale factor is 1.0.

The transfer function is determined to be normalized when its value at
the calper is within 0.01 of 1.0.
The data are determined not to have been scaled if the difference
between calib and scale is < 0.01.
This is because when data are scaled, SAC sets scale to 1.0, but leaves
the calib value olone.
*/

enum NormalizationStatus { Normalized, UnNormalized, Unknown } NormStatus;
enum ScalingStatus { Scaled, UnScaled, ScaleUnknown } ScaleStatus;
sac *s;

double DataMultiplier = 1.0;
DBlist tree; /* used in search for calib, calper */
double calib, calper; /* Used to determine whether data are scaled and
whether response is normalized. */
double EffectiveCalib;
s = sacget_current();
/* If a wfdisc file is loaded, then get the Instrument Information from
there. If not, set to a default value
*/
tree = smGetDefaultTree();

if (!GetCalibCalperFromWfdisc(tree, s->h->nwfid, &calib, &calper)) {
calib = CALIB_UNDEF;
calper = CALPER_UNDEF;
}

/*
Attribute calper
Units ( "Seconds" )
Null ( "-1.0" )
Description ( "nominal calibration period" )
Detail {
This gives the period for which calib, ncalib and calratio
are valid.
}

Attribute calib
Units ( "Nanometers/digital count" )
Null ( "0.0" )
Description ( "nominal calibration" )
Detail {
This is the conversion factor that maps digital data to
displacement, velocity, or acceleration, depending on the
value of segtype or rsptype. The factor holds true at the
oscillation period specified by the attribute calper. A
positive value means ground motion (velocity,
acceleration) increasing in the component direction (up,
north, east) is indicated by increasing counts. A
negative value means the opposite. Calib generally
reflects the best calibration information available at the
time of recording, but refinement may be given in sensor
reflecting a subsequent recalibration of the instrument.
See calratio.
};

Inserted/Modified from passcal/database/data/schemas/css3.0
Your location may vary.

*/
/* Now determine whether response has been normalized or not...
Calper should only be set if a CSS file was read in
*/
if (calper <= 0) {
NormStatus = Unknown;
} else {
NormStatus =
IsNormalized(calper, nfreq, delfrq, xre,
xim) ? Normalized : UnNormalized;
}

/* Now determine whether scale has been applied...
Calib should only be set if a CSS file was read in
*/
if (EQUAL(calib, CALIB_UNDEF)) {
ScaleStatus = ScaleUnknown;
} else {
if (EQUAL(calib, s->h->scale)) {
ScaleStatus = UnScaled;
} else {
ScaleStatus = Scaled;
}
}

EffectiveCalib = calib;
if (NormStatus == UnNormalized && ScaleStatus == Scaled) {
DataMultiplier = 1.0 / EffectiveCalib;
} else if (NormStatus == Normalized && ScaleStatus == UnScaled) {
DataMultiplier = EffectiveCalib;
}
if (!EQUAL(DataMultiplier, 1.0)) {
printf("\tWaveform multiplied by %f after deconvolution.\n",
DataMultiplier);
}

return DataMultiplier;
}

/* -------------------------------------------------------------------- */

void /*FUNCTION*/
transfer(dat, npts, delta, fpfrom, ipfrom, kpfrom, kpfrom_s, fpto, ipto, kpto,
kpto_s, f, iprew, sre, sim, nfft, xre, xim, water_level, nfreq, nerr)
float dat[];
int npts;
double delta;
float fpfrom[];
int ipfrom[];
char kpfrom[MAXKP][MCPFN + 1];
int kpfrom_s;
float fpto[];
int ipto[];
char kpto[MAXKP][MCPFN + 1];
int kpto_s;
double f[];
int *iprew;
double sre[], sim[];
int nfft;
double xre[], xim[];
double water_level;
int nfreq, *nerr;
{
char errmsg[131];
int i;
float a[21];
float nmScale = 1; /* for EVALRESP: scales meters to nanometers. */
double delfrq, denr;
float h_max, threshold;
sac *s;

double DataMultiplier = 1.0;

/*=================================================================
* PURPOSE: To apply an instrument transfer function to a data set.
*=================================================================
* INPUT ARGUMENTS:
*=================================================================
* OUTPUT ARGUMENTS:
*=================================================================
* MODULE/LEVEL: ICM/4
*=================================================================
* SUBROUTINES CALLED:
*=================================================================
* MODIFICATION HISTORY:
* 100325: Cleaned up ranges in do loops
* 871023: Major restructuring of instrument parameter storage.
* 870317: Now passing in work arrays rather than local ones.
* 870219: Added prewhiten/dewhiten for flatter spectrum
* and changed order of spectral operations
* 861105: Original version from K. Nakanishi's TRANSFER program.
*==================================================================
* DOCUMENTED/REVIEWED: 870317
*================================================================ */
s = sacget_current();
/* PROCEDURE: */
delfrq = 1.0e0 / ((double) (nfft) * (double) (delta));

/* - Prewhiten the data if requested. */
if (*iprew > 0) {
errmsg[0] = '\0';
prewit(dat, npts, (float) delta, iprew, a, NULL, errmsg);
if (errmsg[0])
fprintf(stdout, "%s \n", errmsg);
}

/* - Deconvolve seismometer 'FROM' diretion (xre, xim) */

setTransferDirection(FROM);
dseis(nfreq, delfrq, sre, sim, fpfrom, ipfrom, kpfrom, kpfrom_s, &nmScale,
nerr);
if (*nerr != 0)
goto L_8888;

printf(" Station (%s), Channel (%s)\n", s->h->kstnm, s->h->kcmpnm);

/* Attempt to determine whether data have been scaled or not and if
transfer function is normalized or not. Based on that information,
the data may need to be scaled. This function must be called immediately
after getting the deconvolution transfer function before the real and
imaginary arrays have been modified. */

DataMultiplier = GetNormalizationFactor(nfreq, delfrq, sre, sim);

/* compute transfer function
Input: (xre,xim) 'FROM' Transfer direction
Output: (sre, sim) 1.0/ 'FROM' Transfer direction
Set to zero if less than FLT_MIN */

h_max = 0.0f;
for (i = 1; i < nfreq; i++) {
denr = (powi(sre[i], 2) + powi(sim[i], 2));
if (denr > h_max)
h_max = denr;
}
threshold = h_max * pow(10.0, -water_level/20.0);
//printf("water_level = %20.12e\n", water_level);
//printf("h_max = %20.12e, threshold=%20.12e, factor=%20.12e\n", h_max, threshold, pow(10.0, -water_level/20.0));
for (i = 1; i < nfreq; i++) {
denr = (powi(sre[i], 2) + powi(sim[i], 2));
//if (denr <= FLT_MIN) {
if (denr <= threshold) {
sre[i] = 0.0e0;
sim[i] = 0.0e0;
} else {
denr = 1.0e0 / denr;
sre[i] = sre[i] * denr;
sim[i] = -sim[i] * denr;
}
}

/* - Determine seismometer transfer function in the 'TO' direction
(xre, xim). */

setTransferDirection(TO);
dseis(nfreq, delfrq, xre, xim, fpto, ipto, kpto, kpto_s, &nmScale, nerr);

if (*nerr != 0)
goto L_8888;

if((is_type(kpto[0], "none") || is_type(kpto[0], "vel") || is_type(kpto[0],"acc"))) {
int to_type = -1;
if(is_type(kpto[0], "none")) {
to_type = 1;
} else if(is_type(kpto[0], "vel")) {
to_type = 2;
} else if(is_type(kpto[0],"acc")) {
to_type = 3;
}
if(is_type(kpfrom[0], "evalresp")) {
// Here we assume the output units are in Nanometers
switch(to_type) {
case 1: printf(" Units: nm (nanometers)\n"); break;
case 2: printf(" Units: nm/s (nanometers / second)\n"); break;
case 3: printf(" Units: nm/s^2 (nanometers / second^2)\n"); break;
}
} else if(is_type(kpfrom[0], "polezero")) {
if(strcasecmp(kpfrom[2], "meters") == 0 ||
strcasecmp(kpfrom[2], "meter") == 0 ||
strcasecmp(kpfrom[2], "m") == 0) {
switch(to_type) {
case 1: printf(" Units: m (meters)\n"); break;
case 2: printf(" Units: m/s (meters / second)\n"); break;
case 3: printf(" Units: m/s^2 (meters / second^2)\n"); break;
}
} else {
if(strlen(kpfrom[2]) > 0 ) {
printf(" Units: Unrecognized input units in polezero file (%s)\n", kpfrom[2]);
} else {
printf(" Units: Input units not defined in polezero file\n");
}
}
}
}

/*
Apply TO and FROM Instrument responses
F(DATA) * TO
DATA = F^-1(----------------)
FROM
DATA is overwritten
*/

ztransfer(dat, npts, delta, sre, sim, xre, xim, nfreq, nfft, delfrq, f);

/* - Copy the transformed data back into the original data array.
double precieions => single precision */

/* nmScale is 1 by default, but if EVALRESP is used, nmScale converts
to nm. DataMultiplier is used to apply or unapply calib */

for (i = 0; i < npts; i++)
dat[i] = sre[i] * nmScale * DataMultiplier;

/* - Undo the effects of prewhitening if necessary. */

if (*iprew > 0) {
errmsg[0] = '\0';
dewit(dat, npts, *iprew, a, errmsg);
if (errmsg[0])
fprintf(stdout, "%s \n", errmsg);
}

L_8888:
return;

} /* end of function */

#include <string.h>

#include "icm.h"
#include "bool.h"
#include "hdr.h"
#include "amf.h"
#include "dfm.h"

#include "EVRESPnames.h"

#include "bot.h"
#include "ucf.h"
#include "msg.h"
#include "cpf.h"
#include "co.h"
#include "dff.h"
#include "defs.h"

ICM_EXTERN

void DisconnectFromOracleTransfer(void);

void /*FUNCTION*/
xtransfer(nerr)
int *nerr;
{
int ldone;
int iprewu, jdfl, maxn, nfft, nfreq, nra, ntused;

double *sre, *sim, *xre, *xim, water_level;
sac *s;
/*=======================================================================
* PURPOSE: To parse and execute the action command transfer.
* This command deconvolves a seismogram into ground displace-
* ment, then convolves the displacementgram with the response
* of another instrument.
*=======================================================================
* OUTPUT ARGUMENTS:
* nerr: Error flag. Set to 0 if no error occurred. [i]
*=======================================================================
* MODULE/LEVEL: icm/2
*=======================================================================
* GLOBAL INPUT:
* mach: mlarge
* icm: kinstr, ninstr
* dfm: ndfl
*=======================================================================
* GLOBAL OUTPUT:
* icm: fpfrom, lfpfrom, ipfrom, lipfrom, kpfrom, lkpfrom,
* fpto, lfpto, ipto, lipto, kpto, lkpto,
* freq, lfreql, lprew, iprew
* hdr: depmin, depmax, depmen
* mem: sacmem
*=======================================================================
* SUBROUTINES CALLED:
* saclib: lcmore, cfmt, cresp
*=======================================================================
* MODIFICATION HISTORY:
* 871023: Major restructuring of instrument parameter storage.
* 870316: Minor modifications before putting into SAC.
* 861105: Original version in TESTSHELL.
*=======================================================================
* DOCUMENTED/REVIEWED: 871023
*======================================================================= */
/* PROCEDURE: */
*nerr = 0;
xim = NULL;
xre = NULL;
sim = NULL;
sre = NULL;
/* PARSING PHASE: */

/* - Loop on each token in command: */

ldone = FALSE;
clearEVRESPstrucs();

water_level = 600.;
while (lcmore(nerr)) {

/* -- "FROM ...": set 'FROM' instrument information. */
if (lckey("FROM#$", 7)) {
setTransferDirection(FROM);
getins((char *) kmicm.kinstr, 9, cmicm.ninstr, &ldone, cmicm.fpfrom,
cmicm.lfpfrom, cmicm.ipfrom, cmicm.lipfrom,
(char *) kmicm.kpfrom, MCPFN + 1, cmicm.lkpfrom, nerr);
if (*nerr != 0)
goto L_8888;
if (ldone)
break;
}

/* -- "TO ...": set 'TO' instrument information. */
else if (lckey("TO#$", 5)) {
setTransferDirection(TO);
getins((char *) kmicm.kinstr, 9, cmicm.ninstr, &ldone, cmicm.fpto,
cmicm.lfpto, cmicm.ipto, cmicm.lipto, (char *) kmicm.kpto,
MCPFN + 1, cmicm.lkpto, nerr);
if (*nerr != 0)
goto L_8888;
if (ldone)
break;
}

/* -- "FREQLIMITS v1 v2 v3 v4": set frequency limits. */
else if (lkra("FREQ#LIMITS$", 13, 4, 4, cmicm.freq, &nra)) {
if (nra == 4)
cmicm.lfreql = TRUE;
}

/* -- "PREWHITEN ON|OFF|n": set prewhitener order. */
else if (lklogi("PREW#HITEN$", 12, &cmicm.lprew, &cmicm.iprew)) { /* do nothing */
}

/* -- "water_level": set prewhitener order. */
else if (lkreal("WATER#LEVEL$", 13, &water_level)) { /* do nothing */
}

/* -- Bad syntax. */
else {
cfmt("Illegal Option:", 17);
cresp();
}
}
if (water_level < 0.0)
water_level = 0.0;


/* - The above loop is over when one of two conditions has been met:
* (1) An error in parsing has occurred. In this case nerr is > 0 .
* (2) All the tokens in the command have been successfully parsed. */

if (*nerr != 0)
goto L_8888;

/* CHECKING PHASE: */

/* - Test for completeness of input parameters */

/* -- Assume 'NONE' as instrument type if one was omitted. */
if (!cmicm.lkpfrom[0]) {
fstrncpy(kmicm.kpfrom[0], MCPFN, "NONE", 4);
cmicm.lkpfrom[0] = TRUE;
}
if (!cmicm.lkpto[0]) {
fstrncpy(kmicm.kpto[0], MCPFN, "NONE", 4);
cmicm.lkpto[0] = TRUE;
}

/* -- Check the "from" and "to" instrument parameters. */
ckinst(cmicm.fpfrom, cmicm.lfpfrom, cmicm.ipfrom, cmicm.lipfrom,
(char *) kmicm.kpfrom, MCPFN + 1, cmicm.lkpfrom, nerr);
if (*nerr != 0)
goto L_8888;
ckinst(cmicm.fpto, cmicm.lfpto, cmicm.ipto, cmicm.lipto,
(char *) kmicm.kpto, MCPFN + 1, cmicm.lkpto, nerr);
if (*nerr != 0)
goto L_8888;

/* -- Check the frequency limits. */
if ((cmicm.freq[0] >= cmicm.freq[1]) || (cmicm.freq[1] >= cmicm.freq[2]) ||
(cmicm.freq[2] >= cmicm.freq[3])) {
setmsg("WARNING", 2111);
outmsg();
clrmsg();
cmicm.freq[0] = -2.;
cmicm.freq[1] = -1.;
cmicm.freq[2] = 1.e5;
cmicm.freq[3] = 1.e6;
}

/* - Test for a non-null data file list. */

vflist(nerr);

/* - Make sure each file is an evenly spaced time series file. */

vfeven(nerr);
if (*nerr != 0)
goto L_8888;

/* - EXECUTION PHASE: */

/* - Set up prewhitening order. */

if (cmicm.lprew) {
iprewu = cmicm.iprew;
} else {
iprewu = 0;
}

/* Determine the maximum number of data points in the files in DFL. */

vfmax(&maxn, &ntused);

/* - Set up scratch space for TRANSFER subroutine.
* It needs four DOUBLE PRECISION arrays,two for the ffts that must be
* sized to the next power of two above the maximum number of points(nfft)
* and two that are half that size plus 1(nfreq.) */

nfft = next2(maxn);
nfreq = nfft / 2 + 1;
sre = (double *) malloc(sizeof(double) * nfft);
sim = (double *) malloc(sizeof(double) * nfft);
xre = (double *) malloc(sizeof(double) * nfft);
xim = (double *) malloc(sizeof(double) * nfft);
if (!sre || !sim || !xre || !xim) {
goto L_8888;
}

/* - Perform the requested function on each file in DFL. */

for (jdfl = 1; jdfl <= saclen(); jdfl++) {
if (!(s = sacget(jdfl - 1, TRUE, nerr))) {
goto L_8888;
}
//getfil( jdfl, TRUE, &nlen, &ndxy, &ndxx, nerr );

/* -- Call the specific subroutine to work on this file.
* (This is usually a call to a subroutine.) */

transfer(s->y, s->h->npts, DT(s), cmicm.fpfrom, cmicm.ipfrom,
kmicm.kpfrom, MCPFN + 1, cmicm.fpto, cmicm.ipto, kmicm.kpto,
MCPFN + 1, cmicm.freq, &iprewu, sre, sim, nfft, xre, xim, water_level,
nfreq, nerr);

if (*nerr != 0)
goto L_8888;

/* -- Update any header fields that may have changed. */
extrma(s->y, 1, s->h->npts, &s->h->depmin, &s->h->depmax,
&s->h->depmen);

if (kmicm.kpto[0][0]) {
modcase(TRUE, kmicm.kpto[0], strlen(kmicm.kpto[0]), kmicm.kpto[0]);
if (!strncmp(kmicm.kpto[0], "DIS", 3) ||
!strncmp(kmicm.kpto[0], "NONE", 4))
s->h->idep = IDISP;
else if (!strncmp(kmicm.kpto[0], "VEL", 3))
s->h->idep = IVEL;
else if (!strncmp(kmicm.kpto[0], "ACC", 3))
s->h->idep = IACC;
}

if (*nerr != 0)
goto L_8888;

}

/* - Calculate and set new range of dependent variable. */

setrng();

L_8888:
/* - Release scratch space. */
FREE(sre);
FREE(sim);
FREE(xre);
FREE(xim);

DisconnectFromOracleTransfer();
return;

} /* end of function */

Attachments
21:54:42 v.ad6b513c