#include #include #include /* memset() */ #include #include /* fstat() */ #include /* fstat() */ #ifndef COMPLEX_DEF #include "complex.h" #endif #ifndef RAM_UTIL_DEF #include "ram-util.h" #endif AEL * read_ael_output (FILE *fp, int *num_ael) { AEL *ael; char *linebuf; int num_arrays, num_elements; int i, j, k; double tx, ty, tz; double ox, oy, oz; linebuf = (char *) calloc (BUFFER_LENGTH, sizeof(char)); if (linebuf == NULL) { perror ("calloc linebuf in read_ael_output"); return(NULL); } /* Read in first comment line */ if (fgets(linebuf, BUFFER_LENGTH, fp) == NULL) { perror ("fgets linebuf in read_ael_output()"); fprintf (stderr, "error reading comment line from AEL output file\n"); return(NULL); } fprintf (stderr, "Using AEL information: %s\n", linebuf); /* Read in number of arrays */ memset(linebuf, 0, BUFFER_LENGTH); if (fgets(linebuf, BUFFER_LENGTH, fp) == NULL) { perror ("fgets linebuf in read_ael_output()"); fprintf (stderr, "error reading number of arrays line from AEL output file\n"); return(NULL); } if (sscanf (linebuf, "%d", &num_arrays) != 1) { fprintf (stderr, "error parsing num_arrays line: %s", linebuf); fprintf (stderr, "line should be a single integer (like, say, 5)\n"); return(NULL); } /* Read in total number of elements */ memset(linebuf, 0, BUFFER_LENGTH); if (fgets(linebuf, BUFFER_LENGTH, fp) == NULL) { perror ("fgets linebuf in read_ael_output()"); fprintf (stderr, "error reading number of elements line from AEL output file\n"); return(NULL); } if (sscanf (linebuf, "%d", &num_elements) != 1) { fprintf (stderr, "error parsing num_elements line: %s", linebuf); fprintf (stderr, "line should be a single integer (like, say, 150)\n"); return(NULL); } /* Allocate space for positions */ ael = (AEL *) calloc (num_elements, sizeof(AEL)); if (ael == NULL) { perror ("calloc AEL in read_ael_output()"); return(NULL); } /* Read in positions */ k = 0; for (i = 0; i < num_arrays; i++) { for (j = 0; j < (num_elements / num_arrays); j++) { memset(linebuf, 0, BUFFER_LENGTH); if (fgets(linebuf, BUFFER_LENGTH, fp) == NULL) { perror ("fgets linebuf in read_ael_output()"); fprintf (stderr, "error reading ael line from AEL output file\n"); return(NULL); } if (sscanf (linebuf, "%lf %lf %lf", &tx, &ty, &tz) != 3) { fprintf (stderr, "error parsing AEL line: %s", linebuf); return(NULL); } /* Store the offsets only (for now) */ if (j == 0) { ox = tx; oy = ty; oz = tz; } /* Store the results */ ael[k].acoid = i + 1; ael[k].cn = j + 1; ael[k].x = tx - ox; ael[k].y = ty - oy; ael[k].z = oz - tz; /* Increment k */ k++; } /* j loop */ } /* i loop */ /* All done */ *num_ael = num_elements; return(ael); } void correct_ael (AEL *ael, int num_ael, double tide_offset, double **p) { /* AEL values are given WRT a WGS-84 datum. The next step is to * offset for the bottom anchor depth and tide changes. Also, the * receiver depths must be reversed, so zero is at the surface of * the ocean. */ int i; for (i = 0; i < num_ael; i++) { ael[i].z = tide_offset + p[ael[i].acoid - 1][2] - ael[i].z; } } RAM * read_ram_output (FILE *fp) { RAM *ram; /* Reads in the RAM header */ ram = read_ram_output_header (fp); if (ram == NULL) { fprintf (stderr, "error reading RAM output header in read_ram_output()\n"); return(NULL); } /* Read in the RAM data */ if (read_ram_pressure_field (fp, ram) < 0) { fprintf (stderr, "error reading RAM pressure field in read_ram_output()\n"); return(NULL); } return(ram); } RAM * read_ram_output_header (FILE *fp) { RAM *ram; int recsize; int nreadz; int twiddle = 2; int i, j; int lr1 = 0, lr2 = 0; float r; struct stat filestatus; /* * Taken from the "loadramout" matlab script written by Brian Sperry */ /* Allocate return structure */ ram = (RAM *) calloc (1, sizeof(RAM)); if (ram == NULL) { perror ("calloc ram in read_ram_output()"); return(NULL); } /* * Read in the first record */ if (fread(&(ram->lz), sizeof(int), 1, fp) != 1) { perror ("fread ram->lz from RAM output file"); return(NULL); } if (fread(&(ram->lr), sizeof(int), 1, fp) != 1) { perror ("fread ram->lr from RAM output file"); return(NULL); } if (fread(&(ram->f), sizeof(double), 1, fp) != 1) { perror ("fread ram->f from RAM output file"); return(NULL); } /* * Advance to the second record. We're assuming the field is output * here, not TL (and the record length reflects this) */ recsize = 2 * ram->lz * 4 * twiddle; nreadz = 2; if (fseek(fp, 1 * recsize, SEEK_SET) < 0) { perror ("fseek to record 1 in RAM output file"); return(NULL); } /* * Read in the second record */ if (fread(&(ram->zmin), sizeof(double), 1, fp) != 1) { perror ("fread ram->zmin from RAM output file"); return(NULL); } if (fread(&(ram->zmax), sizeof(double), 1, fp) != 1) { perror ("fread ram->zmax from RAM output file"); return(NULL); } if (fread(&(ram->rmin), sizeof(double), 1, fp) != 1) { perror ("fread ram->rmin from RAM output file"); return(NULL); } if (fread(&(ram->rmax), sizeof(double), 1, fp) != 1) { perror ("fread ram->rmax from RAM output file"); return(NULL); } /* * Advance to the next record */ if (fseek(fp, 2 * recsize, SEEK_SET) < 0) { perror ("fseek to record 2 in RAM output file"); return(NULL); } if (ram->lr < 0) { /* Read in the ranges and dr */ if (fread(&(ram->rng0), sizeof(double), 1, fp) != 1) { perror ("fread ram->rng0 from RAM output file"); return(NULL); } if (fread(&(ram->rng1), sizeof(double), 1, fp) != 1) { perror ("fread ram->rng1 from RAM output file"); return(NULL); } if (fread(&(ram->dr), sizeof(double), 1, fp) != 1) { perror ("fread ram->dr from RAM output file"); return(NULL); } /* * Advance to the next record */ if (fseek(fp, 3 * recsize, SEEK_SET) < 0) { perror ("fseek to record 3 in RAM output file"); return(NULL); } if (ram->dr < 0) { /* Altering range step "on the fly" */ if (fread(&(ram->rmax1), sizeof(double), 1, fp) != 1) { perror ("fread ram->rmax1 from RAM output file"); return(NULL); } if (fread(&(ram->dr), sizeof(double), 1, fp) != 1) { perror ("fread ram->dr from RAM output file"); return(NULL); } if (fread(&(ram->dr2), sizeof(double), 1, fp) != 1) { perror ("fread ram->dr2 from RAM output file"); return(NULL); } /* * Advance to the next record */ if (fseek(fp, 4 * recsize, SEEK_SET) < 0) { perror ("fseek to record 3 in RAM output file"); return(NULL); } } } /* * Solve for the number of ranges */ if (ram->lr < 0) { ram->lr = 0; for (r = ram->rmin; r <= ram->rmax; r+=ram->dr) if (r >= ram->rng0 && r <= ram->rng1) (ram->lr)++; if (ram->dr2 != 0) { lr1 = (int) (ram->rmax1 / ram->dr + 0.999); lr2 = (int) ((ram->rmax - ram->rmax1) / ram->dr2 + 0.999); ram->lr = lr1 + lr2; /* Save the grid starting point for future work */ ram->grid_start = lr1 - 1; } ram->dz = (ram->zmax - ram->zmin) / (float) (ram->lz - 1); } /* Allocate space for the z and r indices */ ram->z = (double *) calloc (ram->lz, sizeof(double)); if (ram->z == NULL) { perror ("calloc ram->z in read_ram_output()"); return(NULL); } ram->r = (double *) calloc (ram->lr, sizeof(double)); if (ram->r == NULL) { perror ("calloc ram->r in read_ram_output()"); return(NULL); } /* Fill the z and r indices */ for (i = 0; i < ram->lz; i++) ram->z[i] = ram->zmin + i * ram->dz; if (ram->dr2 == 0) { i = j = 0; while (ram->rmin + j * ram->dr <= ram->rng1) { if (ram->rmin + j * ram->dr >= ram->rng0) ram->r[i++] = ram->rmin + j * ram->dr; j++; } } else { for (i = 0; i < (lr1-1); i++) ram->r[i] = (i+1) * ram->dr; ram->r[lr1-1] = (int) ((ram->rmax1 / ram->dr) + 0.9999) * ram->dr; for (i = lr1; i < (lr1 + lr2 - 1); i++) ram->r[i] = ram->rmax1 + (i+1-lr1) * ram->dr2; ram->r[lr1+lr2-1] = ram->r[lr1-1] + (int) ((ram->rmax - ram->rmax1) / ram->dr2 + 0.999) * ram->dr2; } /* for (i = 0; i < ram->lr; i++) if (ram->rmin + i * ram->dr >= ram->rng0 && ram->rmin + i * ram->dr <= ram->rng1) ram->r[i] = ram->rmin + i * ram->dr; */ /* * Save the current byte offset */ ram->cf_offset = ftell(fp); /* * Allocate a buffer for the RAM data */ ram->cf_buffer = (complex *) calloc (RAM_BUFFER_SIZE, sizeof(complex)); if (ram->cf_buffer == NULL) { perror ("calloc ram->cf_buffer in read_ram_output_header()"); return(NULL); } /* Record the file size */ memset (&filestatus, 0, sizeof(struct stat)); if (fstat(fileno(fp), &filestatus) < 0) { perror ("fstat RAM file in read_ram_output_header()"); return(NULL); } ram->file_size = filestatus.st_size; /* * All done! */ return(ram); } complex * read_recip_pressure_field (int n, FILE **fp, RAM **ram, int rindex, int zindex) { int i; complex *cf; /* Allocate complex pressure vector */ cf = (complex *) calloc(n, sizeof(complex)); if (cf == NULL) { perror ("calloc cf in read_recip_pressure_field()"); return(NULL); } /* Read in the pressures */ for (i = 0; i < n; i++) cf[i] = read_ram_pressure_cell (fp[i], ram[i], ram[i]->grid_start + rindex, zindex); /* All done */ return(cf); } complex read_ram_pressure_cell (FILE *fp, RAM *ram, int rindex, int zindex) { int i, cells_read = 0, cells_to_read = 0; unsigned int requested_offset, requested_cell; /* Calculate the requested offset, in bytes */ requested_cell = ram->lz * rindex + zindex; requested_offset = ram->cf_offset + requested_cell * sizeof(complex); /* Check to see if the requested data falls in the current buffer. If so, return the desired value. */ if (requested_cell >= ram->buffer_start && requested_cell < ram->buffer_end) return(ram->cf_buffer[requested_cell - ram->buffer_start]); /* If the requested offset is larger than the file, return a zero */ if (requested_offset > ram->file_size) return(make_complex (-1.0, 0.0)); /* If the data does not fall within the current buffer, read a new buffer */ /* Seek to the correct place int the file */ if (fseek(fp, requested_offset, SEEK_SET) < 0) { perror ("fseek to requested offset in read_ram_pressure_cell()"); fprintf (stderr, "want: %d (%d) size: %d\n", requested_cell, requested_offset, ram->file_size); fprintf (stderr, "cf_offset: %d lz: %d rindex: %d zindex: %d\n", ram->cf_offset, ram->lz, rindex, zindex); return(make_complex (-1.0, 0.0)); } else { cells_to_read = (requested_offset + RAM_BUFFER_SIZE * sizeof(complex) > ram->file_size) ? (ram->file_size - requested_offset) / sizeof(complex) : RAM_BUFFER_SIZE; /* Fill the buffer */ cells_read = fread(ram->cf_buffer, sizeof(complex), cells_to_read, fp); } if (cells_read < cells_to_read) { fprintf (stderr, "?"); for (i = cells_read; i < RAM_BUFFER_SIZE; i++) ram->cf_buffer[i] = make_complex (0.0, 0.0); } /* Set the buffer extents */ ram->buffer_start = requested_cell; ram->buffer_end = requested_cell + cells_read; /* Return the requested data */ return(ram->cf_buffer[0]); } int read_ram_pressure_field (FILE *fp, RAM *ram) { int i; /* * Seek to the correct place in the file */ if (fseek(fp, ram->cf_offset, SEEK_SET) < 0) { perror ("fseek in read_ram_pressure_field()"); return(-1); } /* Allocate space for the complex pressure field */ ram->cf = (complex **) calloc (ram->lr, sizeof(complex *)); if (ram->cf == NULL) { perror ("calloc ram->cf in read_ram_output()"); return(-1); } for (i = 0; i < ram->lr; i++) { ram->cf[i] = (complex *) calloc (ram->lz, sizeof(complex)); if (ram->cf[i] == NULL) { perror ("calloc ram->cf[i] in read_ram_output()"); return(-1); } } /* * The first complex pressure field written by RAM is for range * dr from the source. We don't want this range, so it * will be skipped. */ if (fread (ram->cf[0], sizeof(complex), ram->lz, fp) != ram->lz) { perror ("fread ram->cf[0] in read_ram_output() [first range]"); return(-1); } /* * Read in complex pressure field */ for (i = 0; i < (ram->lr - 1); i++) { if (fread (ram->cf[i], sizeof(complex), ram->lz, fp) != ram->lz) { perror ("fread ram->cf[i] in read_ram_output()"); fprintf (stderr, "i: %d/%d\n", i, ram->lr); return(-1); } } /* * All done */ return(0); } double ** read_bottom_receiver_positions (FILE *fp) { /* * Read the XYZ bottom HF hydrophone positions in */ char linebuf[BUFFER_LENGTH]; double **p; int acoid; double bx, by, bz; int i; /* * Allocate space for the bottom receiver positions (xyz) */ p = (double **) calloc (NUM_ARRAYS, sizeof(double *)); if (p == NULL) { perror ("calloc p in read_bottom_receiver_positions()"); return(NULL); } for (i = 0; i < NUM_ARRAYS; i++) { p[i] = (double *) calloc (3, sizeof(double)); if (p[i] == NULL) { perror ("calloc p[i] in read_bottom_receiver_positions()"); return(NULL); } } memset (linebuf, 0, BUFFER_LENGTH); if (fgets (linebuf, BUFFER_LENGTH, fp) == NULL) { perror ("fgets"); fprintf (stderr, "error reading in array bottom positions\n"); return(NULL); } while (!feof(fp)) { if (linebuf[0] != '#' && sscanf (linebuf, "%d %lf %lf %lf", &acoid, &bx, &by, &bz) == 4) { if (acoid > 0 && acoid <= 5) { /* store the array positions */ p[acoid-1][0] = bx; p[acoid-1][1] = by; p[acoid-1][2] = bz; } } /* Read in the next line */ memset (linebuf, 0, BUFFER_LENGTH); if (fgets (linebuf, BUFFER_LENGTH, fp) == NULL && !feof(fp)) { perror ("fgets"); fprintf (stderr, "error reading in array bottom positions\n"); return(NULL); } } /* while !feof */ /* All done */ return(p); } int write_ram_output (FILE *fp, RAM *ram, unsigned int writemag) { int i, j; complex_component c; /* Writes output from RAM for perusal with MATLAB. This is intended to be a sanity check... */ fprintf (stderr, "writing %d z and %d r grids\n", ram->lz, ram->lr); if (fwrite (&(ram->lz), sizeof(int), 1, fp) != 1) { perror ("fwrite ram->lz in write_ram_output()"); return(-1); } if (fwrite (&(ram->lr), sizeof(int), 1, fp) != 1) { perror ("fwrite ram->lr in write_ram_output()"); return(-1); } fprintf (stderr, "this run was simulated at %f Hz\n", ram->f); if (fwrite (&(ram->f), sizeof(double), 1, fp) != 1) { perror ("fwrite ram->f in write_ram_output()"); return(-1); } fprintf (stderr, "writing depth and range vectors\n"); if (fwrite(ram->z, sizeof(double), ram->lz, fp) != ram->lz) { perror ("fwrite ram->z in write_ram_output()"); return(-1); } if (fwrite(ram->r, sizeof(double), ram->lr, fp) != ram->lr) { perror ("fwrite ram->r in write_ram_output()"); return(-1); } /* Check to see if phase or magnitude were requested */ if (writemag) { /* Magnitude */ fprintf (stderr, "writing magnitude, columnwise\n"); for (i = 0; i < ram->lr; i++) { for (j = 0; j < ram->lz; j++) { c = complex_magnitude (ram->cf[i][j]); if (fwrite(&c, sizeof(complex_component), 1, fp) != 1) { perror ("fwrite c in write_ram_output()"); fprintf (stderr, "range: %d/%d depth: %d/%d\n", i, ram->lr, j, ram->lz); return(-1); } } /* j loop */ } /* i loop */ } else { /* Phase */ fprintf (stderr, "writing phase, columnwise\n"); for (i = 0; i < ram->lr; i++) { for (j = 0; j < ram->lz; j++) { c = complex_phase (ram->cf[i][j]); if (fwrite(&c, sizeof(complex_component), 1, fp) != 1) { perror ("fwrite c in write_ram_output()"); fprintf (stderr, "range: %d/%d depth: %d/%d\n", i, ram->lr, j, ram->lz); return(-1); } } /* j loop */ } /* i loop */ } /* mag/phase if */ /* All done! */ return(0); } int write_output_file_header (OutputFileHeader ofh, FILE *ofp) { if (fwrite(&(ofh.byteorder), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.byteorder in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.acoid), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.acoid in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.cn), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.cn in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.fs), sizeof(double), 1, ofp) != 1) { perror ("fwrite ofh.fs in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.data_format), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.data_format in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.starting_frame), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_frame in write_output_file_header()"); return(-1); } /* zero_time: time of zero'th frame */ if (fwrite(&(ofh.zero_time.yearday), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.zero_time.yearday in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.zero_time.hour), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.zero_time.hour in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.zero_time.minute), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.zero_time.minute in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.zero_time.second), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.zero_time.second in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.zero_time.usec), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.zero_time.usec in write_output_file_header()"); return(-1); } /* starting_time: time of first frame in file */ if (fwrite(&(ofh.starting_time.yearday), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_time.yearday in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.starting_time.hour), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_time.hour in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.starting_time.minute), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_time.minute in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.starting_time.second), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_time.second in write_output_file_header()"); return(-1); } if (fwrite(&(ofh.starting_time.usec), sizeof(UINT), 1, ofp) != 1) { perror ("fwrite ofh.starting_time.usec in write_output_file_header()"); return(-1); } /* All set! */ return(0); } complex interpolate_complex_pressure (RAM *ram, RAM_COMPANION *rc, double ael_range, double ael_depth) { double range_index, depth_index; int range_index_int, depth_index_int; int i, j; double weight_sum, weight_factor, distance; complex interp_cf, scale; /* Figure out which index has the pressure field output */ if (ael_range < ram->rng0 || ael_range > ram->rng1) { fprintf (stderr, "Requested range is outside the range written to file.\n"); fprintf (stderr, "Requested range: %f meters\n", ael_range); fprintf (stderr, "Range written in file: %f to %f meters\n", ram->rng0, ram->rng1); return(make_complex(-1.0, 0.0)); } range_index = ((ael_range - ram->rng0) / (double) ram->dr); range_index_int = (int) range_index; if (range_index_int + 1 > ram->lr || range_index_int < 0) { fprintf (stderr, "range index out of bounds for interpolation, sorry!\n"); return(make_complex(-1.0, 0.0)); } /* Now solve for the depth index */ if (ael_depth < ram->zmin || ael_depth > ram->zmax) { fprintf (stderr, "Requested interpolation depth is outside the range written\n"); fprintf (stderr, "to file.\n"); fprintf (stderr, "Requested depth: %f meters\n", ael_depth); fprintf (stderr, "Depth range written to file: %f to %f meters\n", ram->zmin, ram->zmax); return(make_complex(-1.0, 0.0)); } depth_index = (ael_depth / (double) ram->dz) - 1; depth_index_int = (int) depth_index; /* We now have 4 points; calculate the weights on them */ weight_sum = 0.0; for (i = range_index_int; i <= range_index_int + 1; i++) { for (j = depth_index_int; j <= depth_index_int + 1; j++) { distance = sqrt ((ram->r[i] - ael_range) * (ram->r[i] - ael_range) + (ram->z[j] - ael_depth) * (ram->z[j] - ael_depth)); weight_factor = 1.0 / distance; weight_sum += weight_factor; } } /* Now come up with a weighted complex pressure function */ interp_cf = make_complex (0.0, 0.0); for (i = range_index_int; i <= range_index_int + 1; i++) { for (j = depth_index_int; j <= depth_index_int + 1; j++) { distance = sqrt ((ram->r[i] - ael_range) * (ram->r[i] - ael_range) + (ram->z[j] - ael_depth) * (ram->z[j] - ael_depth)); weight_factor = 1.0 / distance; interp_cf = complex_add(interp_cf, complex_mult_component (ram->cf[i][j], weight_factor / weight_sum)); } /* depth loop */ } /* range loop */ /* If a scale needs to be placed on the pressure field, do it now */ if (rc->scale_flag == 2) { scale = complex_exp (make_complex(0.0, 2.0 * M_PI * (ram->f / rc->c0) * ael_range)); /* Take into account range spread */ scale = complex_div_component (scale, sqrt(ael_range)); interp_cf = complex_mult (scale, interp_cf); } /* All set */ return(interp_cf); } complex ** regrid_receiver_elements(RAM *ram, RAM_COMPANION *rc) { /* * The purpose of this routine is to interpolate the results of the * propagation run. Here, the source is located at the position of * an acoustic element. The receiver is "marched" through space to * a maximum range. Unfortunately, each acoustic element is at a * slightly different position. This creates an offset for the * entire array as a function of range (since the zero range is * different at each element). * * This routine sets the zero range to be the location of the bottom * element. The offset in range between */ complex **cf_interp; int i, j; double range_grid_offset = 0.0, w1, w2; int rgo_int; /* Allocate space in memory to hold the grid */ cf_interp = (complex **) calloc (ram->lr, sizeof(complex *)); if (cf_interp == NULL) { perror ("calloc cf_interp in regrid_receiver_elements()"); return(NULL); } for (i = 0; i < ram->lr; i++) { cf_interp[i] = (complex *) calloc (ram->lz, sizeof(complex)); if (cf_interp[i] == NULL) { perror ("calloc cf_interp[i] in regrid_receiver_elements()"); return(NULL); } } /* Regrid in range */ /* range_grid_offset = rc->range_offset / ram->dr;*/ rgo_int = (int)range_grid_offset; w1 = range_grid_offset - (double)rgo_int; w2 = 1.0 - w1; if (range_grid_offset > 0.0) { for (i = 0; i < ram->lr; i++) { if ((i - rgo_int) >= 0 && (i - rgo_int) < ram->lr && (i - (rgo_int + 1)) >= 0 && (i - (rgo_int + 1)) < ram->lr) { /* Loop through depths */ for (j = 0; j < ram->lz; j++) { cf_interp[i][j] = complex_add (complex_mult_component (ram->cf[i - rgo_int][j], w2), complex_mult_component (ram->cf[i - (rgo_int + 1)][j], w1)); } /* j loop; depth */ } /* if */ } /* i loop; range */ } else { for (i = 0; i < ram->lr; i++) { if ((i - rgo_int) >= 0 && (i - rgo_int) < ram->lr && (i - (rgo_int - 1)) >= 0 && (i - (rgo_int - 1)) < ram->lr) { /* Loop through depths */ for (j = 0; j < ram->lz; j++) { cf_interp[i][j] = complex_add (complex_mult_component (ram->cf[i - rgo_int][j], w2), complex_mult_component (ram->cf[i - (rgo_int - 1)][j], w1)); } /* j loop; depth */ } /* if */ } /* i loop; range */ } /* else rgo > 0*/ /* Return */ return(cf_interp); } int write_sparse_pressure (FILE *fp, int header_offset_bytes, int nchan, int cn, int lr, int r_index, int lz, int z_index, complex cf) { int total_offset; int record_size; record_size = nchan * sizeof(complex) + sizeof(double) + sizeof(double); /* Calculate total offset */ total_offset = header_offset_bytes + (r_index * lz + z_index) * record_size + sizeof(double) + sizeof(double) + (cn-1) * sizeof(complex); if (fseek(fp, total_offset, SEEK_SET) < 0) { perror ("fseek in write_sparse_pressure()"); fprintf (stderr, "range %d/%d, depth %d/%d, channel %d/%d\n", r_index, lr, z_index, lz, cn, nchan); return(-1); } /* Write complex pressure */ if (fwrite(&(cf), sizeof(complex), 1, fp) < 0) { perror ("fwrite cf in write_sparse_pressure()"); fprintf (stderr, "range %d/%d, depth %d/%d, channel %d/%d\n", r_index, lr, z_index, lz, cn, nchan); return(-1); } /* All done */ return(0); } int dump_pressure_field (FILE *fp, RAM_COMPANION *rc, complex ***cfptr) { int i, j; double trial_range, trial_depth; for (i = 0; i < rc->grid_nr; i++) { for (j = 0; j < rc->grid_nz; j++) { /* Write range and depth pair */ trial_range = rc->grid_r_min + i * rc->grid_delta_r + rc->recip_r; if (fwrite(&trial_range, sizeof(double), 1, fp) != 1) { perror ("fwrite trial_range in dump_pressure_field()"); return(-1); } trial_depth = rc->grid_z_min + j * rc->grid_delta_z; if (fwrite(&trial_depth, sizeof(double), 1, fp) != 1) { perror ("fwrite trial_depth in dump_pressure_field()"); return(-1); } /* Write complex pressure vector */ if (fwrite(cfptr[i][j], sizeof(complex), NUM_ELEMENTS_PER_ARRAY, fp) != NUM_ELEMENTS_PER_ARRAY) { perror ("fwrite cfptr in dump_pressure_field()"); fprintf (stderr, "range cell %d/%d; depth cell %d/%d\n", i, rc->grid_nr, j, rc->grid_nz); } } /* j loop */ } /* i loop */ return(0); } SVP * read_baseline_water_svp (FILE *fp) { /* * Reads the "standard" water SVP from file. * Returns a series of depth, speed pairs */ SVP *s; s = (SVP *) calloc (1, sizeof(SVP)); if (s == NULL) { perror ("calloc s"); return(NULL); } s->z = (double *) calloc (s->nz + 1, sizeof(double)); if (s->z == NULL) { perror ("calloc s->z in read_baseline_water_svp()"); return(NULL); } s->c = (double *) calloc (s->nz + 1, sizeof(double)); if (s->c == NULL) { perror ("calloc s->c in read_baseline_water_svp()"); return(NULL); } while (!feof(fp)) { if (fscanf (fp, "%lf %lf", &(s->z[s->nz]), &(s->c[s->nz])) == 2) { (s->nz)++; s->z = (double *) realloc (s->z, (s->nz + 1) * sizeof(double)); if (s->z == NULL) { perror ("realloc s->z in read_baseline_water_svp()"); return(NULL); } s->c = (double *) realloc (s->c, (s->nz + 1) * sizeof(double)); if (s->c == NULL) { perror ("realloc s->z in read_baseline_water_svp()"); return(NULL); } } } /* while !feof */ /* All done */ return(s); } void write_ram_setup (FILE *fp, Bathymetric *b, SVP *s, double f, double c0, double sx, double sy, double sz, double bx, double by, double range, double range_step, double bearing, unsigned int scale_flag, unsigned int recip_flag, unsigned int project_flag, double recip_r, double recip_base, double grid_min_r) { /* * Writes non-environment related RAM parameters to * file (range, mesh size, etc.) */ int i; double ael_r = 0.0, fudged_range_step = 0.0; /* Write title line */ fprintf (fp, "RAM run, from %f %f %f to %f %f (range: %f bearing: %f), freq %f, flags: s%dr%dp%d\n", sx, sy, sz, bx, by, range, bearing, f, scale_flag, recip_flag, project_flag); /* Write frequency, source and receiver depths. The receiver depth is ignored, so we'll make it the same as the source depth for now */ fprintf (fp, "%f %f %f\n", f, sz, sz); /* Write the maximum range, stepping range and decimation */ if (recip_flag) { fudged_range_step = grid_min_r / range_step; fudged_range_step = grid_min_r / (int) (fudged_range_step+1); } fprintf (fp, "%f %f %d\n", range + BATHYMETRIC_SAMPLE_INTERVAL * 2.0, (recip_flag) ? fudged_range_step : range_step, 1); /* Write maximum depth to calculate to (arbitrarily set to 1000 meters) the depth increment (0.25 meters), depth decimation (1 - for no decimation) depth output (300 meters) */ fprintf (fp, "%f %f %d %f\n", 1000.0, 0.25, 1, 300.0); /* Write reference sound speed (c0 for now), number of terms in rational approximation (8 for now), number of stability constraints (1 for now), and maximum range of stability constraints (0 for now) */ fprintf (fp, "%f %d %d %f\n", c0, 8, 1, 0.0); /* Write no starter field file */ fprintf (fp, "0\n"); /* Write range extraction parameters */ fprintf (fp, "%d %d %d %d %d\n", 1, /* complex pressure field */ scale_flag, /* spreading loss? scale factor? 0: no spreading loss, but cexp 1: spreading loss, with cexp 2: no scaling (no cexp) */ 0, /* no TL line o/p */ 1, /* restrict output to a specific range */ recip_flag); /* recip flag - change steps "on the fly" */ /* Write output range restriction */ if (recip_flag || project_flag) fprintf (fp, "%f %f\n", 1.0, range); else { if (range < BATHYMETRIC_SAMPLE_INTERVAL) fprintf (fp, "%f %f\n", 1.0, range + BATHYMETRIC_SAMPLE_INTERVAL); else fprintf (fp, "%f %f\n", range - BATHYMETRIC_SAMPLE_INTERVAL, range + BATHYMETRIC_SAMPLE_INTERVAL); } if (recip_flag) { /* Write "on the fly" step changes */ ael_r = recip_r - recip_base; fprintf (fp, "%f %f\n", grid_min_r, range_step); } /* Write bathymetric parameters. Bathymetry is calculated from the base of the array. The simulation is done from an actual element, which may not be at the same x,y position as the base of the array. If this is the case, offset the bathymetry by the difference between the base of the array and the array element. */ for (i = 0; i < b->nr; i++) { if (recip_flag) { fprintf (fp, "%f %f\n", (b->bp[i].r - ael_r) < 0.0 ? 0.0 : (b->bp[i].r - ael_r), b->bp[i].z); } else { fprintf (fp, "%f %f\n", b->bp[i].r, b->bp[i].z); } } fprintf (fp, "-1 -1\n"); /* Write SVP block */ for (i = 0; i < b->nr; i++) { /* Write range point */ if (i > 0) { if (recip_flag) { fprintf (fp, "%f\n", (b->bp[i].r - ael_r) < 0.0 ? 0.0 : (b->bp[i].r - ael_r)); } else { fprintf (fp, "%f\n", b->bp[i].r); } } /* Write water SVP */ (void) write_standard_svp(fp, s, b->bp[i].z); /* Write bottom compressional speed */ fprintf (fp, "0 1560.0\n"); fprintf (fp, "%f 1560.0\n", b->bp[i].z); fprintf (fp, "%f 1821.0\n", b->bp[i].z + 89); fprintf (fp, "%f 1862.0\n", b->bp[i].z + 89); fprintf (fp, "%f 2374.0\n", b->bp[i].z + 389); fprintf (fp, "-1 -1\n"); /* Write bottom density information */ fprintf (fp, "0 1.85\n"); fprintf (fp, "%f 1.85\n", b->bp[i].z); fprintf (fp, "%f 1.88\n", b->bp[i].z + 89); fprintf (fp, "%f 2.03\n", b->bp[i].z + 389); fprintf (fp, "-1 -1\n"); /* Write bottom attenuation information */ fprintf (fp, "0 0.18\n"); fprintf (fp, "%f 0.18\n", b->bp[i].z); fprintf (fp, "%f 0.03\n", b->bp[i].z + 89); fprintf (fp, "%f 0.04\n", b->bp[i].z + 389); fprintf (fp, "-1 -1\n"); /* Cycle around */ } /* Cycle through bathymetric points */ /* All done */ return; } void write_standard_svp (FILE *fp, SVP *s, double zb) { /* Writes the standard SVP (as read from file) to the RAM input file */ int i; for (i = 0; i < s->nz; i++) { /* Be sure not to go deeper than the bathymetry */ if (s->z[i] < zb) { fprintf (fp, "%f %f\n", s->z[i], s->c[i]); } } /* Since the bottom is variable, assume the SVP is more or less a straight line to the bottom */ fprintf (fp, "%f %f\n", zb, s->c[s->nz - 1]); /* Write the final -1, -1 */ fprintf (fp, "-1 -1\n"); return; } void tide_compensation (Bathymetric *b, double tide_offset) { int i; /* * Add "tide_offset" to all bathymetric depths */ for (i = 0; i < b->nr; i++) b->bp[i].z += tide_offset; return; } int write_companion_file(FILE *ofp, int acoid, double sx, double sy, double sdepth, double bx, double by, double c0, double range_step, unsigned int grid_flag, unsigned int scale_flag, unsigned int recip_flag, unsigned int project_flag, double grid_r_min, double grid_delta_r, int grid_r_index, int grid_nr, double grid_z_min, double grid_delta_z, int grid_z_index, int grid_nz, double grid_bearing, double recip_base, double recip_r, double cx, double cy, double center_range, int cn) { /* Write ACO ID and channel numbers */ if (fwrite(&acoid, sizeof(int), 1, ofp) != 1) { perror ("fwrite acoid in write_companion_file()"); return(-1); } /* Write source x, y, and depth for this run */ if (fwrite(&sx, sizeof(double), 1, ofp) != 1) { perror ("fwrite sx in write_companion_file()"); return(-1); } if (fwrite(&sy, sizeof(double), 1, ofp) != 1) { perror ("fwrite sy in write_companion_file()"); return(-1); } if (fwrite(&sdepth, sizeof(double), 1, ofp) != 1) { perror ("fwrite sdepth in write_companion_file()"); return(-1); } /* Write receiver x and y for this run */ if (fwrite(&bx, sizeof(double), 1, ofp) != 1) { perror ("fwrite bx in write_companion_file()"); return(-1); } if (fwrite(&by, sizeof(double), 1, ofp) != 1) { perror ("fwrite by in write_companion_file()"); return(-1); } /* Write c0 */ if (fwrite(&c0, sizeof(double), 1, ofp) != 1) { perror ("fwrite c0 in write_companion_file()"); return(-1); } /* Write range step */ if (fwrite(&range_step, sizeof(double), 1, ofp) != 1) { perror ("fwrite range_step in write_companion_file()"); return(-1); } /* Write scale flag */ if (fwrite(&scale_flag, sizeof(unsigned int), 1, ofp) != 1) { perror ("fwrite scale_flag in write_companion_file()"); return(-1); } /* Write grid flag */ if (fwrite(&grid_flag, sizeof(unsigned int), 1, ofp) != 1) { perror ("fwrite grid_flag in write_companion_file()"); return(-1); } /* Write reciprocating flag */ if (fwrite(&recip_flag, sizeof(unsigned int), 1, ofp) != 1) { perror ("fwrite recip_flag in write_companion_file()"); return(-1); } /* Write projection flag */ if (fwrite(&project_flag, sizeof(unsigned int), 1, ofp) != 1) { perror ("fwrite project_flag in write_companion_file()"); return(-1); } if (grid_flag) { /* Write the grid parameters */ if (fwrite(&grid_r_min, sizeof(double), 1, ofp) != 1) { perror ("fwrite grid_r_min in write_companion_file()"); return(-1); } if (fwrite(&grid_delta_r, sizeof(double), 1, ofp) != 1) { perror ("fwrite grid_delta_r in write_companion_file()"); return(-1); } if (fwrite(&grid_r_index, sizeof(int), 1, ofp) != 1) { perror ("fwrite grid_r_index in write_companion_file()"); return(-1); } if (fwrite(&grid_nr, sizeof(int), 1, ofp) != 1) { perror ("fwrite grid_nr in write_companion_file()"); return(-1); } if (fwrite(&grid_z_min, sizeof(double), 1, ofp) != 1) { perror ("fwrite grid_z_min in write_companion_file()"); return(-1); } if (fwrite(&grid_delta_z, sizeof(double), 1, ofp) != 1) { perror ("fwrite grid_delta_z in write_companion_file()"); return(-1); } if (fwrite(&grid_z_index, sizeof(int), 1, ofp) != 1) { perror ("fwrite grid_z_index in write_companion_file()"); return(-1); } if (fwrite(&grid_nz, sizeof(int), 1, ofp) != 1) { perror ("fwrite grid_nz in write_companion_file()"); return(-1); } if (fwrite(&grid_bearing, sizeof(double), 1, ofp) != 1) { perror ("fwrite grid_bearing in write_companion_file()"); return(-1); } } if (recip_flag) { /* Write reciprocating parameters */ if (fwrite (&recip_base, sizeof(double), 1, ofp) != 1) { perror ("fwrite recip_base in write_companion_file()"); return(-1); } if (fwrite (&recip_r, sizeof(double), 1, ofp) != 1) { perror ("fwrite recip_r in write_companion_file()"); return(-1); } if (fwrite (&cx, sizeof(double), 1, ofp) != 1) { perror ("fwrite cx in write_companion_file()"); return(-1); } if (fwrite (&cy, sizeof(double), 1, ofp) != 1) { perror ("fwrite cy in write_companion_file()"); return(-1); } if (fwrite (¢er_range, sizeof(double), 1, ofp) != 1) { perror ("fwrite center_range in write_companion_file()"); return(-1); } if (fwrite (&cn, sizeof(int), 1, ofp) != 1) { perror ("fwrite cn in write_companion_file()"); return(-1); } } /* recip flag */ if (project_flag) { /* Write projection parameters */ if (fwrite (¢er_range, sizeof(double), 1, ofp) != 1) { perror ("fwrite center_range in write_companion_file()"); return(-1); } } /* project flag */ /* All done! */ return(0); } RAM_COMPANION * read_ram_companion_file (FILE *fp) { RAM_COMPANION *rc; rc = (RAM_COMPANION *) calloc (1, sizeof(RAM_COMPANION)); if (rc == NULL) { perror ("calloc rc in read_ram_companion_file()"); return(NULL); } /* Read ACO ID */ if (fread(&(rc->acoid), sizeof(int), 1, fp) != 1) { perror ("fread rc->acoid in read_companion_file()"); return(NULL); } /* Read source x, y, and depth for this run */ if (fread(&(rc->sx), sizeof(double), 1, fp) != 1) { perror ("fread rc->sx in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->sy), sizeof(double), 1, fp) != 1) { perror ("fread rc->sy in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->sdepth), sizeof(double), 1, fp) != 1) { perror ("fread rc->sdepth in read_ram_companion_file()"); return(NULL); } /* Read receiver x and y for this run */ if (fread(&(rc->bx), sizeof(double), 1, fp) != 1) { perror ("fread rc->bx in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->by), sizeof(double), 1, fp) != 1) { perror ("fread rc->by in read_ram_companion_file()"); return(NULL); } /* Read c0 */ if (fread(&(rc->c0), sizeof(double), 1, fp) != 1) { perror ("fread rc->c0 in read_ram_companion_file()"); return(NULL); } /* Read range step */ if (fread(&(rc->range_step), sizeof(double), 1, fp) != 1) { perror ("fread rc->range_step in read_ram_companion_file()"); return(NULL); } /* Read the scale flag */ if (fread(&(rc->scale_flag), sizeof(unsigned int), 1, fp) != 1) { perror ("fread rc->scale_flag in read_ram_companion_file()"); return(NULL); } /* Read grid flag */ if (fread(&(rc->grid_flag), sizeof(unsigned int), 1, fp) != 1) { perror ("fread rc->grid_flag in read_ram_companion_file()"); return(NULL); } /* Read reciprocal flag */ if (fread(&(rc->recip_flag), sizeof(unsigned int), 1, fp) != 1) { perror ("fread rc->recip_flag in read_ram_companion_file()"); return(NULL); } /* Read projection flag */ if (fread(&(rc->project_flag), sizeof(unsigned int), 1, fp) != 1) { perror ("fread rc->recip_flag in read_ram_companion_file()"); return(NULL); } if (rc->grid_flag) { /* Read the grid parameters */ if (fread(&(rc->grid_r_min), sizeof(double), 1, fp) != 1) { perror ("fread rc->grid_r_min in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_delta_r), sizeof(double), 1, fp) != 1) { perror ("fread rc->grid_delta_r in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_r_index), sizeof(int), 1, fp) != 1) { perror ("fread rc->grid_r_index in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_nr), sizeof(int), 1, fp) != 1) { perror ("fread rc->grid_nr in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_z_min), sizeof(double), 1, fp) != 1) { perror ("fread rc->grid_z_min in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_delta_z), sizeof(double), 1, fp) != 1) { perror ("fread rc->grid_delta_z in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_z_index), sizeof(int), 1, fp) != 1) { perror ("fread rc->grid_z_index in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_nz), sizeof(int), 1, fp) != 1) { perror ("fread rc->grid_nz in read_ram_companion_file()"); return(NULL); } if (fread(&(rc->grid_bearing), sizeof(double), 1, fp) != 1) { perror ("fread rc->grid_bearing in read_ram_companion_file()"); return(NULL); } } /* grid flag */ if (rc->recip_flag) { if (fread (&(rc->recip_base), sizeof(double), 1, fp) != 1) { perror ("fread rc->recip_base in read_ram_companion_file()"); return(NULL); } if (fread (&(rc->recip_r), sizeof(double), 1, fp) != 1) { perror ("fread rc->recip_r in read_ram_companion_file()"); return(NULL); } if (fread (&(rc->cx), sizeof(double), 1, fp) != 1) { perror ("fread rc->cx in read_ram_companion_file()"); return(NULL); } if (fread (&(rc->cy), sizeof(double), 1, fp) != 1) { perror ("fread rc->cy in read_ram_companion_file()"); return(NULL); } if (fread (&(rc->center_range), sizeof(double), 1, fp) != 1) { perror ("fread rc->center_range in read_ram_companion_file()"); return(NULL); } if (fread (&(rc->cn), sizeof(int), 1, fp) != 1) { perror ("fread rc->cn in read_ram_companion_file()"); return(NULL); } } /* recip flag */ if (rc->project_flag) { /* Write projection parameter */ if (fread (&(rc->center_range), sizeof(double), 1, fp) != 1) { perror ("fread rc->center_range in read_ram_companion_file()"); return(NULL); } } /* project flag */ /* All done! */ return(rc); } int initialize_regrid_file (FILE *fp, RAM_COMPANION *rc, RAM *ram, unsigned int incore) { int nbytes = 0; int i, j; int nchan = NUM_ELEMENTS_PER_ARRAY; double trial_range, trial_depth; complex *cf; /* Initializes the "regridded" Green's function simulation file */ (void)rewind(fp); /* Write the ACO */ if (fwrite(&(rc->acoid), sizeof(int), 1, fp) != 1) { perror ("fwrite rc->acoid in initialize_regrid_file()"); return(-1); } nbytes += sizeof(int); /* Write the "base" x, y coordinates */ if (fwrite(&(rc->bx), sizeof(double), 1, fp) != 1) { perror ("Fwrite rc->bx in initialize_regrid_file()"); return(-1); } nbytes += sizeof(double); if (fwrite(&(rc->by), sizeof(double), 1, fp) != 1) { perror ("Fwrite rc->by in initialize_regrid_file()"); return(-1); } nbytes += sizeof(double); /* Write the bearing */ if (fwrite(&(rc->grid_bearing), sizeof(double), 1, fp) != 1) { perror ("fwrite rc->bearing in initialize_regrid_file()"); return(-1); } nbytes += sizeof(double); /* Write the number of trial ranges */ if (fwrite(&(rc->grid_nr), sizeof(int), 1, fp) != 1) { perror ("fwrite rc->grid_nr in initialize_regrid_file()"); return(-1); } nbytes += sizeof(int); /* Write the trial ranges */ for (i = 0; i < rc->grid_nr; i++) { trial_range = rc->grid_r_min + i * rc->grid_delta_r; if (rc->recip_flag) trial_range += rc->recip_r; if (fwrite(&trial_range, sizeof(double), 1, fp) != 1) { perror ("fwrite ram trial range in initialize_regrid_file()"); return(-1); } nbytes += sizeof(double); } /* Write the number of trial depths */ if (fwrite(&(rc->grid_nz), sizeof(int), 1, fp) != 1) { perror ("fwrite rc->grid_nz in initialize_regrid_file()"); return(-1); } nbytes += sizeof(int); /* Write the trial depths */ for (i = 0; i < rc->grid_nz; i++) { trial_depth = rc->grid_z_min + i * rc->grid_delta_z; if (fwrite(&trial_depth, sizeof(double), 1, fp) != 1) { perror ("fwrite ram trial depth in initialize_regrid_file()"); return(-1); } nbytes += sizeof(double); } /* Write the number of channels */ if (fwrite(&nchan, sizeof(int), 1, fp) != 1) { perror ("fwrite number of channels in initialize_regrid_file()"); return(-1); } nbytes += sizeof(int); /* Write the size of the header */ nbytes += sizeof(int); if (fwrite(&nbytes, sizeof(int), 1, fp) != 1) { perror ("fwrite nbytes in initialize_regrid_file()"); return(-1); } if (!incore) { /* Allocate a zero complex pressure vector */ cf = (complex *) calloc (nchan, sizeof(complex)); if (cf == NULL) { perror ("calloc cf in initialize_regrid_file()"); return(-1); } /* Cycle through cells */ for (i = 0; i < rc->grid_nr; i++) { for (j = 0; j < rc->grid_nz; j++) { /* Write the range/depth for this cell */ trial_range = rc->grid_r_min + i * rc->grid_delta_r; if (rc->recip_flag) trial_range += rc->recip_r; if (fwrite(&trial_range, sizeof(double), 1, fp) != 1) { perror ("fwrite trial_range in initialize_regrid_file()"); return(-1); } trial_depth = rc->grid_z_min + j * rc->grid_delta_z; if (fwrite(&trial_depth, sizeof(double), 1, fp) != 1) { perror ("fwrite trial_range in initialize_regrid_file()"); return(-1); } /* Write zero complex vector */ if (fwrite(cf, sizeof(complex), nchan, fp) != nchan) { perror ("fwrite cf in initialize_regrid_file()"); fprintf (stderr, "range cell %d/%d; depth cell %d/%d\n", i, ram->lr, j, ram->lz); return(-1); } } } /* All done */ free (cf); } return(nbytes); } complex spline_complex_pressure (RAM *ram, RAM_COMPANION *rc, double ael_range, double ael_depth) { /* Performs a complex spline of Green's function */ complex *splined_depth; double *splined_depth_real, *splined_depth_imag; complex_component *spline_component, *y2; complex_component interp_cf_real, interp_cf_imag; complex interp_cf, scale; int i, j; /* * Spline in the depth direction first */ /* Create a complex depth vector */ splined_depth = (complex *) calloc (ram->lr, sizeof(complex)); if (splined_depth == NULL) { perror ("calloc splined_depth"); return(make_complex(-1.0, 0.0)); } splined_depth_real = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_real == NULL) { perror ("calloc splined_depth_real"); return(make_complex(-1.0, 0.0)); } splined_depth_imag = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_imag == NULL) { perror ("calloc splined_depth_imag"); return(make_complex(-1.0, 0.0)); } /* build up a new "r" vector at the splined depth */ spline_component = (complex_component *) calloc (ram->lz, sizeof(complex_component)); if (spline_component == NULL) { perror ("calloc spline_component"); return(make_complex(-1.0, 0.0)); } for (i = 0; i < ram->lr; i++) { for (j = 0; j < ram->lz; j++) spline_component[j] = complex_real (ram->cf[i][j]); y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on real part, bailing\n"); return(make_complex(-1.0, 0.0)); } splined_depth_real[i] = real_splint (ram->z, spline_component, y2, ram->lz, ael_depth); free(y2); for (j = 0; j < ram->lz; j++) spline_component[j] = complex_imag (ram->cf[i][j]); y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on imag part, bailing\n"); return(make_complex(-1.0, 0.0)); } splined_depth_imag[i] = real_splint (ram->z, spline_component, y2, ram->lz, ael_depth); free(y2); splined_depth[i] = make_complex(splined_depth_real[i], splined_depth_imag[i]); } /* * Now spline in the range direction */ y2 = real_spline (ram->r, splined_depth_real, ram->lr, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on real range, bailing\n"); return(make_complex(-1.0, 0.0)); } interp_cf_real = real_splint (ram->r, splined_depth_real, y2, ram->lr, ael_range); free (y2); y2 = real_spline (ram->r, splined_depth_imag, ram->lr, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on imag range, bailing\n"); return(make_complex(-1.0, 0.0)); } interp_cf_imag = real_splint (ram->r, splined_depth_imag, y2, ram->lr, ael_range); free (y2); interp_cf = make_complex (interp_cf_real, interp_cf_imag); free (spline_component); free (splined_depth_real); free (splined_depth_imag); free (splined_depth); /* If a scale needs to be placed on the pressure field, do it now */ if (rc->scale_flag == 2) { scale = complex_exp (make_complex(0.0, 2.0 * M_PI * (ram->f / rc->c0) * ael_range)); /* Take into account range spread */ scale = complex_div_component (scale, sqrt(ael_range)); interp_cf = complex_mult (scale, interp_cf); } return(interp_cf); } complex * spline_complex_pressure_range_precision (RAM *ram, RAM_COMPANION *rc, double ael_depth) { /* Performs a complex spline of Green's function - but only in the range direction */ complex *splined_depth; double *splined_depth_real, *splined_depth_imag, *y2r, *y2i; complex_component interp_cf_real, interp_cf_imag; complex_component *spline_component, *y2; complex *interp_cf, scale; int i, j; /* Create a complex depth vector */ splined_depth = (complex *) calloc (ram->lr, sizeof(complex)); if (splined_depth == NULL) { perror ("calloc splined_depth in spline_complex_pressure_range_precision()"); return(NULL); } splined_depth_real = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_real == NULL) { perror ("calloc splined_depth_real in spline_complex_pressure_range_precision()"); return(NULL); } splined_depth_imag = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_imag == NULL) { perror ("calloc splined_depth_imag in spline_complex_pressure_range_precision()"); return(NULL); } interp_cf = (complex *) calloc (rc->grid_nr, sizeof(complex)); if (interp_cf == NULL) { perror ("calloc interp_cf in spline_complex_pressure_range_precision()"); return(NULL); } /* build up a new "r" vector at the splined depth */ spline_component = (complex_component *) calloc (ram->lz, sizeof(complex_component)); if (spline_component == NULL) { perror ("calloc spline_component in spline_complex_pressure_range_precision()"); return(NULL); } for (i = 0; i < ram->lr; i++) { for (j = 0; j < ram->lz; j++) spline_component[j] = complex_real (ram->cf[i][j]); y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on real part, bailing\n"); return(NULL); } splined_depth_real[i] = real_splint (ram->z, spline_component, y2, ram->lz, ael_depth); free(y2); for (j = 0; j < ram->lz; j++) spline_component[j] = complex_imag (ram->cf[i][j]); y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99); if (y2 == NULL) { fprintf (stderr, "error performing real spline on imag part, bailing\n"); return(NULL); } splined_depth_imag[i] = real_splint (ram->z, spline_component, y2, ram->lz, ael_depth); free(y2); splined_depth[i] = make_complex(splined_depth_real[i], splined_depth_imag[i]); } /* * Now spline in the range direction */ y2r = real_spline (ram->r, splined_depth_real, ram->lr, 9.0e99, 9.0e99); if (y2r == NULL) { fprintf (stderr, "error performing real spline on real range, bailing\n"); return(NULL); } y2i = real_spline (ram->r, splined_depth_imag, ram->lr, 9.0e99, 9.0e99); if (y2i == NULL) { fprintf (stderr, "error performing real spline on imag range, bailing\n"); return(NULL); } for (j = 0; j < rc->grid_nr; j++) { interp_cf_real = real_splint (ram->r, splined_depth_real, y2r, ram->lr, rc->grid_r_min + j * rc->grid_delta_r); interp_cf_imag = real_splint (ram->r, splined_depth_imag, y2i, ram->lr, rc->grid_r_min + j * rc->grid_delta_r); interp_cf[j] = make_complex (interp_cf_real, interp_cf_imag); if (rc->scale_flag == 2) { scale = complex_exp (make_complex(0.0, 2.0 * M_PI * (ram->f / rc->c0) * (rc->grid_r_min + j * rc->grid_delta_r))); /* Take into account range spread */ scale = complex_div_component (scale, sqrt(rc->grid_r_min + j * rc->grid_delta_r)); interp_cf[j] = complex_mult (scale, interp_cf[j]); } } /* j loop */ free (splined_depth_real); free (splined_depth_imag); free (splined_depth); free (spline_component); free (y2r); free (y2i); /* If a scale needs to be placed on the pressure field, do it now */ return(interp_cf); } complex * spline_complex_pressure_range (RAM *ram, RAM_COMPANION *rc, double ael_depth) { /* Performs a complex spline of Green's function - but only in the range direction */ complex *splined_depth; double *splined_depth_real, *splined_depth_imag, *y2r, *y2i; complex_component interp_cf_real, interp_cf_imag; complex *interp_cf, scale; int i, j, depth_index; double depth_error; /* * Spline in the range direction only */ /* Create a complex depth vector */ splined_depth = (complex *) calloc (ram->lr, sizeof(complex)); if (splined_depth == NULL) { perror ("calloc splined_depth"); return(NULL); } splined_depth_real = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_real == NULL) { perror ("calloc splined_depth_real"); return(NULL); } splined_depth_imag = (double *) calloc (ram->lr, sizeof(double)); if (splined_depth_imag == NULL) { perror ("calloc splined_depth_imag"); return(NULL); } interp_cf = (complex *) calloc (rc->grid_nr, sizeof(complex)); if (interp_cf == NULL) { perror ("calloc interp_cf in spline_complex_pressure_range()"); return(NULL); } depth_error = fabs(ael_depth - ram->z[0]); depth_index = 0; for (i = 1; i < ram->lz; i++) { if (fabs(ael_depth - ram->z[i]) < depth_error) { depth_index = i; depth_error = fabs(ael_depth - ram->z[i]); } } if (depth_error > 0.0001) fprintf (stderr, "depth error: %f\n", depth_error); for (i = 0; i < ram->lr; i++) { splined_depth_real[i] = complex_real(ram->cf[i][depth_index]); splined_depth_imag[i] = complex_imag(ram->cf[i][depth_index]); splined_depth[i] = make_complex(splined_depth_real[i], splined_depth_imag[i]); } /* * Now spline in the range direction */ y2r = real_spline (ram->r, splined_depth_real, ram->lr, 9.0e99, 9.0e99); if (y2r == NULL) { fprintf (stderr, "error performing real spline on real range, bailing\n"); return(NULL); } y2i = real_spline (ram->r, splined_depth_imag, ram->lr, 9.0e99, 9.0e99); if (y2i == NULL) { fprintf (stderr, "error performing real spline on imag range, bailing\n"); return(NULL); } for (j = 0; j < rc->grid_nr; j++) { interp_cf_real = real_splint (ram->r, splined_depth_real, y2r, ram->lr, rc->grid_r_min + j * rc->grid_delta_r); interp_cf_imag = real_splint (ram->r, splined_depth_imag, y2i, ram->lr, rc->grid_r_min + j * rc->grid_delta_r); interp_cf[j] = make_complex (interp_cf_real, interp_cf_imag); if (rc->scale_flag == 2) { scale = complex_exp (make_complex(0.0, 2.0 * M_PI * (ram->f / rc->c0) * (rc->grid_r_min + j * rc->grid_delta_r))); /* Take into account range spread */ scale = complex_div_component (scale, sqrt(rc->grid_r_min + j * rc->grid_delta_r)); interp_cf[j] = complex_mult (scale, interp_cf[j]); } } /* j loop */ free (splined_depth_real); free (splined_depth_imag); free (splined_depth); free (y2r); free (y2i); /* If a scale needs to be placed on the pressure field, do it now */ return(interp_cf); } void project_ael (double cx, double cy, double **p, double sbearing, AEL *ael, int num_ael) { /* * Projects AEL coordinates onto a new bearing */ int i; double ux, uy, ax, ay; /* Create a unit vector in the direction of the new bearing */ ux = cos(sbearing); uy = sin(sbearing); for (i = 0; i < num_ael; i++) { /* * Subtract to get the relative position from the center of the array */ ax = (ael[i].x + p[ael[i].acoid - 1][0]) - cx; ay = (ael[i].y + p[ael[i].acoid - 1][1]) - cy; /* * Perform a dot product to project the AEL coordinate onto the bearing line, * and store it in the AEL structure */ ael[i].r = ax * ux + ay * uy; } /* Done */ return; }