gtkseisviewgl_test.c

segy 显示程序！希望能给正在做这部分朋友提供一部分资料

/* 
 * gtkseisviewgl_test - Test for rendering seismic data in
 * GtkSeisViewGl widget (synthetic 2D seismc data is created on
 * every run for this purpose)
 *
 * Copyright (C) 2005-2006 Vladimir Bashkardin
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more av.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 *
 * Author:  Vladimir Bashkardin  <[e-mail]>
 */

#include <math.h>
#include <gtk/gtk.h>
#include <gtkseisviewgl/gtkseisviewgl.h>
#include <gtkseisviewgl/gseisviewtoolmotion.h>
#include <gtkseisviewgl/gseisviewtoolhlight.h>
#include <gtkseisviewgl/gseisviewtoolscale.h>
#include <gtkseisviewgl/gseisviewtoolmagnify.h>
#include <gtkseisviewgl/gseisviewaxisz.h>
#include <gtkseisviewgl/gseisviewaxistraces.h>

#ifndef GTK_SEIS_VIEW_GL_THREADED
#define GTK_SEIS_VIEW_GL_THREADED
#endif

static gboolean image_progress (GtkSeisViewGl *seis_view, GdkPixbuf *pixbuf,
                                gpointer data) {
    return TRUE;
}

static void image_rendered (GtkSeisViewGl *seis_view, GdkPixbuf *pixbuf,
                            gpointer data) {
    GError *save_error = NULL;

    if (pixbuf) {
        g_printf ("Writing to file...\n");
        gdk_pixbuf_save (pixbuf, "/tmp/gsegyview_image.png", "png",
                         &save_error, NULL);
        g_object_unref (pixbuf);
        if (save_error) {
            g_printf ("Save error: %s\n", save_error->message);
            g_error_free (save_error);
        }
        g_printf ("Done.\n");
    }
}

static gboolean window_close (GtkWidget *widget, GdkEvent* event, gpointer data) {
#ifdef GTK_SEIS_VIEW_GL_THREADED
    gtk_seis_view_gl_stop_rendering_thread (GTK_SEIS_VIEW_GL (data));
#endif
    gtk_main_quit ();
    return FALSE;
}

static gboolean shared_window_close (GtkWidget *widget, GdkEvent* event, gpointer data) {
    return TRUE;
}

static void magnify_viewer_area (GtkWidget *widget, gdouble slow_start_data, gdouble slow_size_data,
                                                    gdouble fast_start_data, gdouble fast_size_data, gpointer data) {
    GtkSeisViewGl *shared_seis_view_gl = GTK_SEIS_VIEW_GL (data);
    gtk_seis_view_gl_set_data_viewport (shared_seis_view_gl, slow_start_data, slow_size_data,
                                                             fast_start_data, fast_size_data);
    gtk_seis_view_gl_redraw (shared_seis_view_gl);
}

static gfloat berlague_impulse (gfloat tau, gfloat t) {
    return (exp (-1000 * (tau - t) * (tau - t)) * sin (2.0 * M_PI * 35.0 * (tau - t) + M_PI / 2));
}

static void calculate_seismogram (gfloat *data, gint trace_num, gint profile_length,
                                  gint start_offset, gint samples_num, gfloat dt,
                                  gint layers_num, gfloat *vp, gfloat *h, gfloat *ro) {
    gint MAX_ANGLE = 75;
    gint offset;
    gint trace_dist; /* distance between traces */
    gfloat *vs; /* array[layers_num] of shear velocities */
    gfloat *pr; /* array[layers_num] of Poisson's ratios */
    gfloat *A; /* array[layers_num - 1] of A in Suey's approximation */
    gfloat *B; /* array[layers_num - 1] of B in Suey's approximation */
    gfloat *ang; /* array[layers_num - 1] of incidence angles */
    gfloat *tm; /* array[layers_num - 1] of one-way times in each layer */
    gfloat *refl; /* array[2 * (layers_num - 1)] of reflection times and amplitudes */
    gfloat tt; /* total one-way time to a reflector */
    gint i, j, k, l; /* simple counters */
    gfloat amp; /* Seismic amplitude */
    gfloat temp;

    trace_dist = profile_length / (gfloat)trace_num;

    vs = (gfloat*)g_malloc (layers_num * sizeof (gfloat));
    pr = (gfloat*)g_malloc (layers_num * sizeof (gfloat));
    for (i = 0; i < layers_num; i++) { /* Calculate Poisson's ratio */
        vs[i] = vp[i] * 0.7;
        temp = (vp[i] * vp[i]) / (vs[i] * vs[i]);
        pr[i] = (temp - 2.0) / (2.0 * temp - 2.0);
    }

    A = (gfloat*)g_malloc ((layers_num - 1) * sizeof (gfloat));
    B = (gfloat*)g_malloc ((layers_num - 1) * sizeof (gfloat));
    for (i = 1; i < layers_num; i++) { /* Calculate Shuey's approximation members */
        A[i - 1] = ((vp[i] - vp[i - 1]) / ((vp[i] + vp[i - 1]) / 2.0) + (ro[i] - ro[i - 1]) / ((ro[i] + ro[i - 1]) / 2.0)) / 2.0;
        temp = ((vp[i] - vp[i - 1]) / ((vp[i] + vp[i - 1]) / 2.0)) /
               ((vp[i] - vp[i - 1]) / ((vp[i] + vp[i - 1]) / 2.0) + (ro[i] - ro[i - 1]) / ((ro[i] + ro[i - 1]) / 2.0));
        temp = temp - 2.0 * (1.0 + temp) * ((1 - (pr[i] + pr[i - 1])) / (1 - (pr[i] + pr[i - 1]) / 2.0));
        B[i - 1] = A[i - 1] * temp + (pr[i] - pr[i - 1]) / (1 - (pr[i] + pr[i - 1]) / 2.0) * (1 - (pr[i] + pr[i - 1]) / 2.0);
    }

    ang = (gfloat*)g_malloc ((layers_num - 1) * sizeof (gfloat));
    tm = (gfloat*)g_malloc ((layers_num - 1) * sizeof (gfloat));
    refl = (gfloat*)g_malloc (2 * (layers_num - 1) * sizeof (gfloat));
    for (i = 0; i < trace_num; i++) { /* Produce trace by trace */
        offset = start_offset + i * trace_dist;
        for (j = 0; j < (layers_num - 1); j++) { /* Cycle through horizons */
            amp = 1.0;
            for (k = 0; k <= MAX_ANGLE * 100; k++) { /* Adjust reflection angle */
                ang[j] = M_PI * (k / 100.0) / 180.0;
                tm[j] = sqrt (h[j] * h[j] + h[j] * tan (ang[j]) * h[j] * tan (ang[j])) / vp[j];
                tt = tm[j];
                temp = fabs (offset) / 2.0 - h[j] * tan (ang[j]);
                for (l = j - 1; l >= 0; l--) { /* Calculate incidence angles and corresponding offsets */
                    if (temp < -(trace_dist / 2.0))
                        break;
                    ang[l] = asin (vp[l] * sin (ang[l + 1]) / vp[l + 1]);
                    temp -= h[l] * tan (ang[l]);
                    tm[l] = sqrt (h[l] * h[l] + h[l] * tan (ang[l]) * h[l] * tan (ang[l])) / vp[l];
                    tt += tm[l];
                }
                if (fabs (temp) < (trace_dist / 2.0))
                    break;
            }
            for (k = 0; k <= j; k++) { /* Do dynamic features */
                amp /= 2.0 * vp[k] * tm[k]; /* Spreading */
                if (k != j) /* Transmission */
                    amp *= (1 - A[k] * A[k]);
                else /* Reflection */
                    amp *= (A[k] + B[k] * sin (ang[k]) * sin (ang[k]));
            }
            tt *= 2.0;
            refl[j * 2] = tt;
            refl[j * 2  + 1] = 1.0/*amp*/;
        }
        for (j = 0; j < samples_num; j++) {
            data[i * samples_num + j] = 0.0;
            for (k = 0; k < (layers_num - 1); k++)
                data[i * samples_num + j] += refl[k * 2 + 1] * berlague_impulse (j * dt, refl[k * 2]);
        }
    }

    g_free (refl);
    g_free (tm);
    g_free (ang);
    g_free (B);
    g_free (A);
    g_free (pr);
    g_free (vs);
}

int main (int argc, char *argv[]) {
    gint trace_num = 96;
    gint start_offset = 0/*-1175*/;
    gint profile_length = 2400;
    gint samples_num = 2048;
    gfloat dt = 0.002;