/*************************************************************** TIMSLIP This program was begun in July 1988 and used initially for determining the calibration factors of the WEPOCS II Moana Wave ADCP. Since the DAS then was based on the HP 9816 and recorded time only to the nearest minute, there were unknown discrepancies between the ADCP time and true time. This time offset causes glitches in the calculation of absolute velocity, primarily during ship accelerations (including turns). Therefore, before trying to calibrate the ADCP, we attempt to determine the time offset at each major acceleration (ship going on or off station) when there was GPS data. These offsets are then applied before the calibration factors are calculated. The time offset is determined by trying to minimize some measure of the acceleration glitch. These measures can include the variance of the absolute reference layer velocity on the ensemble time grid (every 5 minutes) or on the fix grid (about every minute for raw GPS), or it can be something like the difference between the absolute reference layer velocity during the acceleration and the average of its values before and after the acceleration. As of Sun 07-31-1988, this program is still in an experimental state. Parameters and options are set in the source code rather than via a control file. It is expected that components of this program may have more general usefulness. This may require rewriting some of the routines for greater generality or modularity. Wed 08-10-1988 changed to control file input. Thu 01-05-1989 switched to the read_ref from refabs, which filters out BADFLOAT values. 89/01/10 Added the ability to work with course changes as well as speed changes. This requires a change in the control file structure. Fri 12-15-1989 Added ability to specify a fixed time shift in seconds in place of "yes" or "no" in response to question "use_shifted_times"; Added output columns for delta-u and delta-v, to be used in determining the effect, if any, of Schuler oscillations, presumed to be proportional to the North acceleration. A corresponding change was made in ADCPCAL.M CONTROL FILE STRUCTURE: fix_file_type: { simple or HIG } fix_file: reference_file: output_file: year_base= { year base for decimal day calculations } min_n_fixes= { minimum number of fixes for calculation of time shift to proceed. } { The jump in velocity will be between indices (n_refs-1)/2 and (n_refs-1)/2 +1. For example, it will be between 6 and 7 if n_refs is 13. Note that all these indices start from 0, C-style, not from 1. However, the lowest index specified must be >= 1, since the time of the start of ensemble 1 is (approximately) the time recorded with ensemble 0.} n_refs= { number of ensembles used in calculation of time shift } i_ref_l0= { first and last indices of ensembles used } i_ref_l1= { in calibration calculation } i_ref_r0= { "l" is before jump, "r" is after } i_ref_r1= { minimum is 1, not zero } up_thresh= { m/s, for jump detection } down_thresh= turn_speed= { m/s, for turn detection } turn_thresh= { degrees } dtmax= 360 { max seconds between ensembles, to screen gaps } tolerance= 5.e-5 { for BRENT, in days } grid: { fix or ensemble; calculate reference layer velocities averaged between fix times, or over ensemble intervals} use_shifted_times? { yes or no or fixed time shift in seconds; use the best-fit time shift in the calibration } */ #include "common.h" #include "dbext.h" #include "use_db.h" #include "ioserv.h" #include "misc.h" #include "vstat.h" #include "cal.h" /* FIX_TYPE, VELOCITY_TYPE, fix_to_velocity(), ... */ #define TIMSLIP #define TSHIFTMAX 10.0/MINUTES_PER_DAY #define NRMAX 120 #define NFMAX 120 #if PROTOTYPE_ALLOWED static int interp_fix(FIX_TYPE *f); int calc_ship(VELOCITY_TYPE *refs,int i_ref_0,int i_ref_1,VELOCITY_TYPE *ship); double abs_ref_grid(VELOCITY_TYPE *refs); int average_velocity(VELOCITY_TYPE *varray, VELOCITY_TYPE *vavg, int n, double t0, double dtmax); double abs_ref_grid(VELOCITY_TYPE *refs); double shifted_abs_ref_grid(double t); double shifted_vrefabs(double t); double spike_size(int i0, int i1, int i2, int i3, VELOCITY_TYPE *refs); int is_jump(VELOCITY_TYPE refs[], int n_refs, double up_thresh, double down_thresh); int is_turn(VELOCITY_TYPE refs[], int n_refs, double turn_speed, double turn_thresh); int get_ref(FILE *fp_ref, VELOCITY_TYPE refs[], int n_refs); int get_fixes(double t0, double t1); int calc_ddop(VELOCITY_TYPE *refs); int calc_dship(VELOCITY_TYPE *refs); int adcp_calib(VELOCITY_TYPE *dop, VELOCITY_TYPE *ship, double *amp, double *phase); double vrefabs(VELOCITY_TYPE *arefs, FIX_TYPE *afixes, VELOCITY_TYPE *arefabs, int an_ref, int an_fix, int *an_refabs); #else static int interp_fix(); int calc_ship(); double abs_ref_grid(); int average_velocity(); double abs_ref_grid(); double shifted_abs_ref_grid(); double shifted_vrefabs(); double spike_size(); int is_jump(); int is_turn(); int get_ref(); int get_fixes(); int calc_ddop(); int calc_dship(); int adcp_calib(); double vrefabs(); #endif /*-------------------------------------------------------------- Global variables for this module: */ static int (*read_fix_ptr)(); FILE *fp_ref, *fp_fix, *fp_out; static int startup = 1; VELOCITY_TYPE refs[NRMAX]; static VELOCITY_TYPE absrefs[NRMAX]; /* on ref grid */ static VELOCITY_TYPE shifted_refs[NRMAX]; int n_refs; static FIX_TYPE fixes[NFMAX]; static int n_fixes; static VELOCITY_TYPE refabs[NFMAX]; /* on fix grid */ static int n_refabs; int year_base; double up_thresh, down_thresh; double turn_speed, turn_thresh; /* ref indices of left, right ranges for calib: */ static int i_ref_l0, i_ref_l1, i_ref_r0, i_ref_r1; static VELOCITY_TYPE dship, ddop; static double dtmax; /* reject longer intervals */ static UNISTAT_TYPE refstat; /* for calculating variance */ /*-------------------------------------------------------------- Functions: */ /* interp_fix: input with f->t set to desired time; at return, f->x and f->y are interpolated. returns: 1 on error (interpolation time out of array range) 0 otherwise This function is very specific to this application, and quite dumb. It always does a simple linear search from the beginning of the array. */ #if PROTOTYPE_ALLOWED static int interp_fix(FIX_TYPE *f) #else static int interp_fix(f) FIX_TYPE *f; #endif { extern FIX_TYPE fixes[]; extern int n_fixes; int i1 = 1; int i0; double df, dt; while ( i1 < n_fixes && fixes[i1].t < f->t ) i1++; i0 = i1 - 1; if (i1 == n_fixes || fixes[i0].t > f->t) return 1; /* error */ dt = fixes[i1].t - fixes[i0].t; df = fixes[i1].x - fixes[i0].x; f->x = fixes[i0].x + df/dt * (f->t - fixes[i0].t); df = fixes[i1].y - fixes[i0].y; f->y = fixes[i0].y + df/dt * (f->t - fixes[i0].t); return 0; } /* calc_ship: calculates the average ship velocity over the reference velocity ensembles i_ref_0 to i_ref_1. The former must be greater than zero, since the time of the first ensemble is the recorded time of the previous ensemble. RETURNS: 0 on success 1 otherwise */ #if PROTOTYPE_ALLOWED int calc_ship(VELOCITY_TYPE *refs, int i_ref_0, int i_ref_1, VELOCITY_TYPE *ship) #else int calc_ship(refs, i_ref_0, i_ref_1, ship) VELOCITY_TYPE *refs; int i_ref_0, i_ref_1; VELOCITY_TYPE *ship; #endif { FIX_TYPE pos0, pos1; /* interpolated positions */ pos0.t = refs[i_ref_0-1].t; if(interp_fix(&pos0)) return 1; pos1.t = refs[i_ref_1].t; if(interp_fix(&pos1)) return 1; if(fix_to_velocity(&pos0, &pos1, ship)) return 1; else return 0; } /* calc_ship() */ /*-------------------------------------------------------------- FUNCTION: average_velocity PURPOSE: given an array of ADCP ensemble velocities, calculate the average over n ensembles. t0 is the start time of the first ensemble, usually taken as the recorded time of the previous ensemble. dtmax is usually set to a little more than the ensemble duration; if any interval between ensembles exceeds it, the function terminates with an error. RETURNS: 0 if successful 1 otherwise ----------------------------------------------------------------*/ #if PROTOTYPE_ALLOWED int average_velocity(VELOCITY_TYPE *varray, VELOCITY_TYPE *vavg, int n, double t0, double dtmax) #else int average_velocity(varray, vavg, n, t0, dtmax) VELOCITY_TYPE *varray, *vavg; int n; double t0, dtmax; #endif { int i; double dt; vavg->u = 0; vavg->v = 0; dt = varray[0].t - t0; for (i=0; i dtmax) return(1); /* error */ vavg->u += varray[i].u * dt; vavg->v += varray[i].v * dt; } dt = varray[n-1].t - t0; vavg->t = 0.5 * (varray[n-1].t + t0); vavg->u = vavg->u / dt; vavg->v = vavg->v / dt; return 0; } /* average_velocity() */ /*-------------------------------------------------------------- FUNCTION: abs_ref_grid Calculates the absolute reference layer velocity for each ADCP ensemble, assuming the ship velocity is constant between fixes. Note that this does not make much sense unless the interval between fixes is smaller than or equal to the ensemble interval. PARAMETERS: refs = the array of VELOCITY_TYPE containing the reference layer velocities RETURNS: variance of the reference layer velocities. --------------------------------------------------------------*/ #if PROTOTYPE_ALLOWED double abs_ref_grid(VELOCITY_TYPE *refs) #else double abs_ref_grid(refs) VELOCITY_TYPE *refs; #endif { VELOCITY_TYPE ship; double variance; int i; /* initialize the statistics routines */ allocate_unistat(&refstat, 2, "abs ref layer", U_VARIANCE); zero_unistat(&refstat); for (i=1; i1, the ADCP is underestimating the magnitude of the flow; correct by multiplying ADCP velocity by amplitude. phase (output) = pointer to misalignment angle in degrees; correct by rotating the measured velocity by this angle. The full correction to a measured velocity expressed as a complex number, u+iv, is U+iV = (u+iv) * amp * exp(i*phase) RETURN: 1 error 0 success ----------------------------------------------------------------*/ #if PROTOTYPE_ALLOWED int adcp_calib(VELOCITY_TYPE *dop, VELOCITY_TYPE *ship, double *amp, double *phase) #else int adcp_calib(dop, ship, amp, phase) VELOCITY_TYPE *dop, *ship; double *amp, *phase; #endif { double uuvv, duv, u2v2d, a_denom; duv = (dop->v * ship->u) - (dop->u * ship->v); uuvv = (dop->u * ship->u) + (dop->v * ship->v); u2v2d = sqr(dop->u) + sqr(dop->v); if (uuvv == 0.0) return 1; *phase = -atan( duv / uuvv ); if ((a_denom = (cos(*phase) * u2v2d)) == 0.0) return 1; *amp = uuvv / a_denom; /* The sign on *amp is opposite to that given by Pollard and Read, because here we are using velocity of ship relative to water, instead of the reverse. The only effect is on the sign of the amplitude. The sign of the phase is reversed so that it agrees with the convention of Joyce (1988) and with the way Ramon's original rotate program worked: The angle is measured counterclockwise from the gyro axis to the transducer axis. */ *phase *= (180.0 / M_PI); return 0; } #if PROTOTYPE_ALLOWED int get_fixes(double t0, double t1) #else int get_fixes(t0, t1) double t0, t1; #endif { int finished; static int n; /* The pair of fixes bracketing the last time of the previous call become the first pair of fixes to test for bracketing of the first time of the present call. */ if (startup) { finished = (*read_fix_ptr)(fp_fix, &(fixes[0])); finished |= (*read_fix_ptr)(fp_fix, &(fixes[1])); n = 2; startup = 0; } /* Make sure the first fix is before t0. If not, t0 cannot be bracketed, so there is no point in looking through more fixes. The program should look for a later velocity jump instead. */ if (fixes[0].t > t0) return (0); /* Look first in the existing array of fixes for the pair bracketing the first time, t0. */ while ( (n>2) && !( (fixes[0].t < t0) && (fixes[1].t >= t0) ) ) { memcpy(&(fixes[0]), &(fixes[1]), (n-1)*sizeof(FIX_TYPE)); n--; } /* If the bracketing pair is not in memory, get it from the fix file. */ while ( !( (fixes[0].t <= t0) && (fixes[1].t > t0) )) { fixes[0] = fixes[1]; finished = (*read_fix_ptr)(fp_fix, &(fixes[1])); if (finished) return (EOF); /* end of file exit */ } /* Now the first time is bracketed by fixes */ while ( !(fixes[n-1].t > t1) && (n < NFMAX)) { finished = (*read_fix_ptr)(fp_fix, &(fixes[n])); if (finished) return (EOF); /* end of file exit */ n++; } return (n); } /* get_fixes() */ /*============================================================*/ #if PROTOTYPE_ALLOWED void do_it(FILE *fp_cnt) #else void do_it(fp_cnt) FILE *fp_cnt; #endif { int finished = 0; int i; int jump_dir; /* 1 is up, -1 is down in speed */ double ref_var; /* reference layer variance */ double ref_var0, ref_var1; /* with - and + TSHIFTMAX */ double ref_varc; /* with const_tshift_d */ double min_ref_var; double brent(); /* function that minimizes a function of one variable */ static double tol = 5.e-5, /* about 5 seconds */ xmin; /* parameters for brent */ double amp, phase; /* for calibration */ static FILE_NAME_TYPE ref_filename, fix_filename, out_filename; static VELOCITY_TYPE *selected_refs; static int use_shift; static char shift_spec[20]; static double const_tshift_d; static int const_tshift = 0; static int min_n_fixes; static int i_fix_file; /* index into array of function pointers */ static int i_grid; /* index into array of function pointers */ /*--------------------------------------------------------- Order of elements must be coordinated between the following two arrays, since the second names the index into the first. */ static double (*shift_var_func_ptr[])()= { shifted_vrefabs, shifted_abs_ref_grid }; static NAME_LIST_ENTRY_TYPE grid_names[]= { {"fix", 0}, {"ensemble", 1}, {NULL, 0} }; static NAME_LIST_ENTRY_TYPE yes_no[]= { {"no", 0}, {"yes", 1}, {NULL, 0} }; /*------------------------------------------------------*/ extern NAME_LIST_ENTRY_TYPE fix_file_types[]; extern int (*read_fix_ptr_array[])(); static NAME_LIST_TYPE name_list[]= { fix_file_types, grid_names, yes_no }; static PARAMETER_LIST_ENTRY_TYPE parameters[] = { {" fix_file_type:", " %s", &i_fix_file, TYPE_TRANS}, {" fix_file:", " %79s", fix_filename, TYPE_STRING}, {" reference_file:", " %79s", ref_filename, TYPE_STRING}, {" output_file:", " %79s", out_filename, TYPE_STRING}, {" year_base=", " %d", &year_base, TYPE_INT}, {" min_n_fixes=", " %d", &min_n_fixes, TYPE_INT}, {" n_refs=", " %d", &n_refs, TYPE_INT}, {" i_ref_l0=", " %d", &i_ref_l0, TYPE_INT}, {" i_ref_l1=", " %d", &i_ref_l1, TYPE_INT}, {" i_ref_r0=", " %d", &i_ref_r0, TYPE_INT}, {" i_ref_r1=", " %d", &i_ref_r1, TYPE_INT}, {" up_thresh=", " %lf", &up_thresh, TYPE_DOUBLE}, {" down_thresh=", " %lf", &down_thresh, TYPE_DOUBLE}, {" turn_speed=", " %lf", &turn_speed, TYPE_DOUBLE}, {" turn_thresh=", " %lf", &turn_thresh, TYPE_DOUBLE}, {" dtmax=", " %lf", &dtmax, TYPE_DOUBLE}, {" tolerance=", " %lf", &tol, TYPE_DOUBLE}, {" grid:", " %s", &i_grid, TYPE_TRANS + 256}, {" use_shifted_times?"," %s", shift_spec, TYPE_STRING}, { NULL, NULL, NULL, 0 } }; /* parameters[] */ if (get_parameters(fp_cnt, parameters, name_list)) exit(-1); if (i_ref_l0 < 1 || i_ref_l1 < i_ref_l0 || i_ref_r0 < i_ref_l1 || i_ref_r1 < i_ref_r0) { printf("\nError in i_ref indices; minimum value is 1, they must increase.\n"); exit(-1); } if ( (use_shift = get_code(yes_no, shift_spec)) == BADINT) /* not yes or no */ { check_error( (const_tshift = atoi(shift_spec)) == 0, "\nuse_shifted_times? , , or ; is shift in sec.\n"); use_shift = 1; const_tshift_d = (double)const_tshift / (double)SECONDS_PER_DAY; } selected_refs = (use_shift ? shifted_refs : refs); /* If the ensemble grid is used, n_refabs will always be the same as n_refs; otherwise, n_refabs will be set for each calibration point by shifted_vrefabs(). */ n_refabs = n_refs; read_fix_ptr = read_fix_ptr_array[i_fix_file]; fp_ref = check_fopen(ref_filename, "r"); fp_fix = check_fopen(fix_filename, "r"); fp_out = check_fopen(out_filename, "w"); /* initialize by filling the refs array and reading two fixes */ for (i=0; i min_n_fixes) { ref_var = (*shift_var_func_ptr[i_grid])(0.0); ref_var1 = (*shift_var_func_ptr[i_grid])(TSHIFTMAX); ref_var0 = (*shift_var_func_ptr[i_grid])(-TSHIFTMAX); #if 0 fprintf(fp_out, "variance of %d absolutes: %f\n", n_refabs); fprintf(fp_out, " + TSHIFTMAX: %f, - TSHIFTMAX: %f\n", ref_var1 = (*shift_var_func_ptr[i_grid])(TSHIFTMAX), ref_var0 = (*shift_var_func_ptr[i_grid])(-TSHIFTMAX)); fprintf(fp_out, " + 1 min: %f, - 1 min: %f\n", (*shift_var_func_ptr[i_grid])( 1.0/MINUTES_PER_DAY), (*shift_var_func_ptr[i_grid])(-1.0/MINUTES_PER_DAY)); #endif if (ref_var < ref_var1 && ref_var < ref_var0) { min_ref_var = brent(-TSHIFTMAX, 0.0, TSHIFTMAX, shift_var_func_ptr[i_grid], tol, &xmin); if (const_tshift != 0) { ref_varc = (*shift_var_func_ptr[i_grid])(const_tshift_d); } if (calc_ddop(selected_refs)) printf("\nerror in calc_ddop"); else if (calc_dship(selected_refs)) printf("\nerror in calc_dship"); else if (adcp_calib(&ddop, &dship, &, &phase)) printf("\nerror in adcp_calib"); else /* everything is ok, so print results on a single line */ { fprintf(fp_out, " %2d %7.3f %2d %3d %6.2f %6.2f %7.3f %7.3f", jump_dir, (refs[(n_refs-1)/2].t), n_refabs, n_fixes, dship.u, dship.v, ref_var, min_ref_var); fprintf(fp_out, " %4.0f %7.3f %7.3f\n", xmin*SECONDS_PER_DAY, amp, phase); } } } } } fclose(fp_ref); fclose(fp_fix); fclose(fp_out); } /* do_it() */