Skip to content

Commit

Permalink
Velocity models (#6)
Browse files Browse the repository at this point in the history
* add velocity model object and allocate_velocity_model utility function

* .cc file to .cpp

* new velocity testing file

* set up separate files for velocity and libseis

* simple velocity model function for generating quick velocity model with just a starting velocity and gradient
  • Loading branch information
dpgraham4401 authored May 5, 2024
1 parent aff5fb5 commit c8bb196
Show file tree
Hide file tree
Showing 9 changed files with 236 additions and 153 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -34,4 +34,6 @@ include/*
lib/*
bin/*
test/test_runner
Testing
test/Temporary
*.bin
9 changes: 7 additions & 2 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,14 @@ set(SRC_DIR "src")

### create libseis library
include_directories(src)
add_library(${LIB_NAME} SHARED ${SRC_DIR}/libseis.c)
add_library(${LIB_NAME} SHARED ${SRC_DIR}/libseis.c ${SRC_DIR}/velocity.c)


set_target_properties(
${LIB_NAME} PROPERTIES
VERSION "v${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}"
PUBLIC_HEADER ${SRC_DIR}/libseis.h
PUBLIC_HEADER ${SRC_DIR}/velocity.h
)

install(
Expand All @@ -41,7 +43,10 @@ set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
FetchContent_MakeAvailable(googletest)

enable_testing()
add_executable(libseis_test test/libseis_test.cc)
add_executable(libseis_test test/libseis_test.cpp)
add_executable(velocity_test test/velocity_test.cpp)
target_link_libraries(libseis_test ${LIB_NAME} GTest::gtest_main)
target_link_libraries(velocity_test ${LIB_NAME} GTest::gtest_main)
include(GoogleTest)
gtest_discover_tests(libseis_test)
gtest_discover_tests(velocity_test)
108 changes: 55 additions & 53 deletions src/libseis.c
Original file line number Diff line number Diff line change
@@ -1,75 +1,77 @@
#include "libseis.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

double *read_double(char path[], int nt, int nx) {
FILE *fptr;
double *data;
data = (double *) malloc(nt * nx * sizeof(double));
fptr = fopen(path, "rb");
if (fptr == NULL) {
printf("unable to open file at %s\n", path);
exit(EXIT_FAILURE);
}
for (int i = 0; i < nt * nx; i++) {
fread(&data[i], sizeof(double), 1, fptr);
}
fclose(fptr);
return data;
FILE *fptr;
double *data;
data = (double *)malloc(nt * nx * sizeof(double));
fptr = fopen(path, "rb");
if (fptr == NULL) {
printf("unable to open file at %s\n", path);
exit(EXIT_FAILURE);
}
for (int i = 0; i < nt * nx; i++) {
fread(&data[i], sizeof(double), 1, fptr);
}
fclose(fptr);
return data;
}

void write_double(char path[], double *data, int nt, int nx) {
FILE *fptr;
fptr = fopen(path, "wb");
if (fptr == NULL) {
printf("unable to open file at %s\n", path);
exit(EXIT_FAILURE);
}
for (int i = 0; i < nx * nt; i++) {
fwrite(&data[i], sizeof(double), 1, fptr);
}
fclose(fptr);
FILE *fptr;
fptr = fopen(path, "wb");
if (fptr == NULL) {
printf("unable to open file at %s\n", path);
exit(EXIT_FAILURE);
}
for (int i = 0; i < nx * nt; i++) {
fwrite(&data[i], sizeof(double), 1, fptr);
}
fclose(fptr);
}

/**
* Creates a new gather with the same data as the input gather, but with each sample
* multiplied by a power of its corresponding time value.
* Creates a new gather with the same data as the input gather, but with each
* sample multiplied by a power of its corresponding time value.
*
* @param gather The input gather to be multiplied.
* @param power The power to raise the time values to.
* @return A new gather with the same dimensions as the input gather, but with modified data.
* @return A new gather with the same dimensions as the input gather, but with
* modified data.
*/
Gather *gain_gather(Gather *gather, double power) {
Gather *new_gather = malloc(sizeof(Gather));
if (new_gather == NULL) {
fprintf(stderr, "Error: input gather is NULL\n");
}
memcpy(new_gather, gather, sizeof(Gather));
new_gather->data = malloc(gather->nt * gather->nx * sizeof(double));
if (new_gather->data == NULL) {
fprintf(stderr, "Error: input gather data is NULL\n");
}
memcpy(new_gather->data, gather->data, gather->nt * gather->nx * sizeof(double));
for (int trace_index = 0; trace_index < gather->nx; trace_index++) {
for (int sample_index = 0; sample_index < gather->nt; sample_index++) {
double time = (double) sample_index * new_gather->dt;
double t_pow = pow(time, power);
int data_index = (trace_index * new_gather->nt) + sample_index;
new_gather->data[data_index] = gather->data[data_index] * t_pow;
}
Gather *new_gather = malloc(sizeof(Gather));
if (new_gather == NULL) {
fprintf(stderr, "Error: input gather is NULL\n");
}
memcpy(new_gather, gather, sizeof(Gather));
new_gather->data = malloc(gather->nt * gather->nx * sizeof(double));
if (new_gather->data == NULL) {
fprintf(stderr, "Error: input gather data is NULL\n");
}
memcpy(new_gather->data, gather->data,
gather->nt * gather->nx * sizeof(double));
for (int trace_index = 0; trace_index < gather->nx; trace_index++) {
for (int sample_index = 0; sample_index < gather->nt; sample_index++) {
double time = (double)sample_index * new_gather->dt;
double t_pow = pow(time, power);
int data_index = (trace_index * new_gather->nt) + sample_index;
new_gather->data[data_index] = gather->data[data_index] * t_pow;
}
return new_gather;
}
return new_gather;
}

void display_gather(Gather *gather) {
printf("--------------------\n");
printf("Gather: id %d\n", gather->id);
printf("Sample Rate: %f seconds\n", gather->dt);
printf("Sampling Frequency: %f Hz\n", 1 / gather->dt);
printf("Nyquist Frequency: %f Hz\n", 1 / (gather->dt * 2));
printf("Samples per trace: %d\n", gather->nt);
printf("Number of traces: %d\n", gather->nx);
printf("--------------------\n");
printf("--------------------\n");
printf("Gather: id %d\n", gather->id);
printf("Sample Rate: %f seconds\n", gather->dt);
printf("Sampling Frequency: %f Hz\n", 1 / gather->dt);
printf("Nyquist Frequency: %f Hz\n", 1 / (gather->dt * 2));
printf("Samples per trace: %d\n", gather->nt);
printf("Number of traces: %d\n", gather->nx);
printf("--------------------\n");
}
43 changes: 24 additions & 19 deletions src/libseis.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,37 +3,40 @@
// David Graham - 2023-15-01

typedef enum {
SINGLE,
DOUBLE,
SINGLE,
DOUBLE,
} PRECISION;

/**
* @brief A collection of seismic traces and attributes
*
* @details A gather holds encapsulates all the information about seismic gathers,
* including the number of traces, the time sampling, number of time samples, and the data itself.
* It's intended to be generic enough to be used by various seismic gather types, such as CMP, or shot.
* @details A gather holds encapsulates all the information about seismic
* gathers, including the number of traces, the time sampling, number of time
* samples, and the data itself. It's intended to be generic enough to be used
* by various seismic gather types, such as CMP, or shot.
*
* @param
* id: (int) The gather's unique identifier
* @param
* nt: (int) number times samples in each trace, assumed homogeneous number of samples in each trace
* nt: (int) number times samples in each trace, assumed homogeneous number of
* samples in each trace
* @param
* nx: (int) number traces
* @param
* dt: (double) the time sampling of the traces, assumed homogeneous
* @param
* data: (pointer) pointer to the array of samples
* @param
* precision: (enum PRECISION) the precision of the samples (e.g., double, double)
* precision: (enum PRECISION) the precision of the samples (e.g., double,
* double)
*/
typedef struct Gather {
int id;
int nt;
int nx;
double dt;
double *data;
PRECISION precision;
int id;
int nt;
int nx;
double dt;
double *data;
PRECISION precision;
} Gather;

double *read_double(char path[], int nt, int nx);
Expand All @@ -44,19 +47,21 @@ void write_double(char path[], double *data, int nt, int nx);
* gain_gather
*
* @details
* Applies a time-variant scaling to a series of seismic traces in the same gather.
* It can be used to compensate for exponential decay (constant Q) in amplitude with time
* See "Tpow: an estimator of seismic amplitude decay" by Fowler & Claerbout
* nt' = nt * (dt * t^pow)
* Applies a time-variant scaling to a series of seismic traces in the same
* gather. It can be used to compensate for exponential decay (constant Q) in
* amplitude with time See "Tpow: an estimator of seismic amplitude decay" by
* Fowler & Claerbout nt' = nt * (dt * t^pow)
*
* @param
* gather: A pointer to a shot, CMP, or other seismic gather with all necessary fields.
* gather: A pointer to a shot, CMP, or other seismic gather with all necessary
* fields.
* @param
* power: The power (double) applied to the seismic trace as a function of time.
*
* @returns
* A pointer to a new gather with the supplied Tpow applied.
* All gather attributes, besides the seismic data, are copied from the original gather
* All gather attributes, besides the seismic data, are copied from the original
* gather
*
*/
Gather *gain_gather(Gather *gather, double power);
Expand Down
28 changes: 28 additions & 0 deletions src/velocity.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
#include "velocity.h"

#include <stdlib.h>

VelocityModel *allocate_velocity_model(const int rows, const int cols) {
VelocityModel *model = malloc(sizeof(VelocityModel));
model->rows = rows;
model->cols = cols;

model->data = (double **)malloc(rows * sizeof(double *));
for (int i = 0; i < rows; i++) {
model->data[i] = (double *)malloc(cols * sizeof(double));
}

return model;
}

VelocityModel *simple_velocity_model(VelocityModel *model,
const double gradient,
const double startingVelocity) {

for (int i = 0; i < model->rows; i++) {
for (int j = 0; j < model->cols; j++) {
model->data[i][j] = startingVelocity + gradient * i;
}
};
return model;
}
19 changes: 19 additions & 0 deletions src/velocity.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
//
// Created by dg on 5/4/24.
//

#ifndef VELOCITY_H
#define VELOCITY_H

typedef struct {
int rows;
int cols;
double **data;
} VelocityModel;

VelocityModel *allocate_velocity_model(int rows, int cols);

VelocityModel *simple_velocity_model(VelocityModel *model, double gradient,
double startingVelocity);

#endif // VELOCITY_H
79 changes: 0 additions & 79 deletions test/libseis_test.cc

This file was deleted.

Loading

0 comments on commit c8bb196

Please sign in to comment.