Convert to double precision (#5)

* convert all current work to double as default * add initial test for gain function * add display_gather function for fun
dpgraham4401 · Dec 4, 2022 · fe6dbd3 · fe6dbd3
1 parent ae49abf
commit fe6dbd3
Show file tree

Hide file tree

Showing 3 changed files with 92 additions and 72 deletions.
diff --git a/src/libseis.c b/src/libseis.c
@@ -4,49 +4,60 @@
 #include <string.h>
 #include <math.h>
 
-float *read_float(char path[], int nt, int nx) {
+double *read_double(char path[], int nt, int nx) {
     FILE *fptr;
-    float *data;
-    data = (float *) malloc(nt * nx * sizeof(float));
+    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(float), 1, fptr);
+        fread(&data[i], sizeof(double), 1, fptr);
     }
     fclose(fptr);
     return data;
 }
 
-void write_float(char path[], float *data, int nt, int nx) {
+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(float), 1, fptr);
+        fwrite(&data[i], sizeof(double), 1, fptr);
     }
     fclose(fptr);
 }
 
-Gather *gain_gather(Gather *gather, float pow) {
+Gather *gain_gather(Gather *gather, double power) {
     Gather *new_gather = malloc(sizeof(struct Gather));
     new_gather->id = gather->id;
     new_gather->nt = gather->nt;
     new_gather->nx = gather->nx;
     new_gather->dt = gather->dt;
-    new_gather->data = malloc(gather->nt * gather->nt * sizeof(float));
+    new_gather->data = malloc(gather->nt * gather->nt * sizeof(double));
     for (int trace_index = 0; trace_index < gather->nx; trace_index++) {
         for (int sample_index = 0; sample_index < gather->nt; sample_index++) {
-            float time = (float) sample_index * new_gather->dt;
-            float t_pow = powf(time, pow);
+            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;
 }
+
+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");
+}
diff --git a/src/libseis.h b/src/libseis.h
@@ -3,8 +3,8 @@
 // David Graham - 2023-15-01
 
 typedef enum {
-  SINGLE,
-  DOUBLE,
+    SINGLE,
+    DOUBLE,
 } PRECISION;
 
 /**
@@ -21,23 +21,24 @@ typedef enum {
  * @param
  * nx: (int) number traces
  * @param
- * dt: (float) the time sampling of the traces, assumed homogeneous
+ * 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., float, double)
+ * precision: (enum PRECISION) the precision of the samples (e.g., double, double)
  */
 typedef struct Gather {
-  int id;
-  int nt; /// a is an int
-  int nx;
-  float dt;
-  float *data;
-  PRECISION precision;
+    int id;
+    int nt;
+    int nx;
+    double dt;
+    double *data;
+    PRECISION precision;
 } Gather;
 
-float *read_float(char path[], int nt, int nx);
-void write_float(char path[], float *data, int nt, int nx);
+double *read_double(char path[], int nt, int nx);
+
+void write_double(char path[], double *data, int nt, int nx);
 
 /**
  * gain_gather
@@ -51,11 +52,13 @@ void write_float(char path[], float *data, int nt, int nx);
  * @param
  * gather: A pointer to a shot, CMP, or other seismic gather with all necessary fields.
  * @param
- * pow: The power (float) applied to the seismic trace as a function of time.
+ * 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
  *
  */
-Gather *gain_gather(Gather *gather, float pow);
+Gather *gain_gather(Gather *gather, double power);
+
+void display_gather(Gather *gather);
diff --git a/test/libseis_test.cc b/test/libseis_test.cc
@@ -1,73 +1,79 @@
 #include "gtest/gtest.h"
 
 extern "C" {
+#include <math.h>
 #include "libseis.h"
 }
 
 int const nt_test = 10;
 int const nx_test = 1;
 
-float float_trace1[nt_test] = {
-    0.000,
-    2.000,
-    4.000,
-    8.000,
-    1.000,
-    2.000,
-    3.000,
-    4.000,
-    5.000,
-    6.000,
+double double_trace1[nt_test] = {
+        0.000,
+        2.000,
+        4.000,
+        8.000,
+        1.000,
+        2.000,
+        3.000,
+        4.000,
+        5.000,
+        6.000,
 };
 
-float cmp_array[nt_test * 2] = {
-    0.000,
-    2.000,
-    4.000,
-    8.000,
-    1.000,
-    2.000,
-    3.000,
-    4.000,
-    5.000,
-    6.000,
-    0.000,
-    0.000,
-    1.000,
-    0.020,
-    0.301,
-    -0.101,
-    0.035,
-    0.011,
-    -0.006,
-    0.009,
+double cmp_array[nt_test * 2] = {
+        0.000,
+        2.000,
+        4.000,
+        8.000,
+        1.000,
+        2.000,
+        3.000,
+        4.000,
+        5.000,
+        6.000,
+        0.000,
+        0.000,
+        1.000,
+        0.020,
+        0.301,
+        -0.101,
+        0.035,
+        0.011,
+        -0.006,
+        0.009,
 };
 
 struct Gather cmp = {
-    1,
-    nt_test,
-    2,
-    1,
-    cmp_array,
-    DOUBLE
+        1,
+        nt_test,
+        2,
+        1,
+        cmp_array,
+        DOUBLE
 };
 
-TEST(HelloTest, ExampleTest) {
-    EXPECT_EQ(7 * 6, 42);
-}
-
 TEST(LibseisIO, ReadWriteFloat) {
     char tmp_file[] = "/tmp/data.bin";
-    write_float(tmp_file, float_trace1, nt_test, nx_test);
-    float *data = read_float(tmp_file, nt_test, nx_test);
+    write_double(tmp_file, double_trace1, nt_test, nx_test);
+    double *data = read_double(tmp_file, nt_test, nx_test);
     for (int i = 0; i < nt_test; i++) {
-        EXPECT_EQ(data[i], float_trace1[i]);
+        EXPECT_EQ(data[i], double_trace1[i]);
     }
 }
 
 TEST(GainCmp, GainsCmpGather) {
-    float t_pow = 2;
+    double t_pow = 2;
     Gather *foo = gain_gather(&cmp, t_pow);
-    printf("%d\n", foo->nt);
-    EXPECT_EQ(0, 0);
+
+    printf("Check that the t^power function operates correctly");
+    for (int i = 0; i < foo->nt; i++) {
+        printf("%f * (%f^%f) = %f\n", cmp.data[i], foo->dt * i, t_pow, foo->data[i]);
+        double t_powered = cmp.data[i] * pow(foo->dt * i, t_pow);
+        EXPECT_EQ(foo->data[i], t_powered);
+    }
+}
+
+TEST(Gather, GatherDisplay) {
+    display_gather(&cmp);
 }