segy_util.c

这是matlab在地球物理数据处理方面的源码

/********************************************************************
* Definition of SEGY headers and access routines to disk SEGY files *
*                                                                   *
* FILE *segy_open_read (char *,                                     *
*                       struct asciiheader *,                       *
*                       struct binaryheader *)                      *
* --- opens file (1st parameter = name),                            *
*     returns ASCII and binary headers and file pointer             *
*     filepointer = NULL for non-SEGY files (already closed)        *
*                                                                   *
* FILE *segy_open_write(char *,                                     *
*                       struct asciiheader *,                       *
*                       struct binaryheader *)                      *
* --- opens file (1st parameter = name),                            *
*     writes ASCII and binary headers, returns file pointer         *
*     check filepointer for value NULL (not open) !                 *
*                                                                   *
* long segy_read (FILE *,                                           *
*                 struct traceheader *,                             *
*                 float data[],                                     *
*                 struct binaryheader *);                           *
* --- reads one trace from file according to binary header          *
*     (used values : format, lproc)                                 *
*     returns trace header, and data converted to float             *
*     return value = sample count (as in traceheader)               *
*                                                                   *
* long segy_write(FILE *,                                           *
*                 struct traceheader *,                             *
*                 float data[])                                     *
* *** only 3 parameters ***                                         *
* --- writes one trace to file according to trace header            *
*     (used value : lproc)                                          *
*     traceheader followed by data converted from float to REAL*4   *
*     return value = sample count (as in traceheader)               *
*                                                                   *
* FILE *image_open_write()                                          *
* long image_write()                                                *
* --- called like the segy... functions. Produce disk image output  *
*     suitable for ProMAX.                                          *
*                                                                   *
*-------------------------------------------------------------------*
*                                                                   *
* ALL parameters are to be submitted by POINTERS !                  *
*                                                                   *
* (C) SGSoft 1992 Version 1.11 920211                               *
*                         1.21 930218 SEGY disk image               *
*                                                                   *
********************************************************************/

#include <stdio.h>
#include <math.h>
#include "segy.h"


/* convert IEEE float to IBM FORTRAN REAL*4 */

void float2real(float xfloat, char *xreal)
{
 long vorzeichen, exponent, mantisse, i;
 double ymantisse;
 union
 {
   unsigned long l;
   char c[4];
 } yreal;

 if (xfloat != 0.0)
 {
  vorzeichen=(xfloat < 0.0);
  xfloat=fabs((double)xfloat);
  exponent=ceil(log((double)xfloat)/log(16.0)+1e-6);
  if (abs(exponent) < 64)
  {
   ymantisse=xfloat/pow(16.0,(double)exponent);
   mantisse=floor(ymantisse*16777216.0); /* 2^24 */
   yreal.l=mantisse | ((exponent+64) << 24) | (vorzeichen << 31);
  }
  else
   yreal.l=0;
 }
 else
  yreal.l=0;
/* poor man's memcpy */
 for (i=0; i<4; i++)
  xreal[i]=yreal.c[i];
}

/* convert IBM FORTRAN REAL*4 to IEEE float */
float real2float(char *xreal)
{
 union
 {
  unsigned long l;
  char c[4];
 } yreal;
 long vorzeichen, exponent, mantisse, i;

 for (i=0; i<4; i++)
  yreal.c[i]=xreal[i];
 vorzeichen=((yreal.l & 0x80000000) != 0);
 exponent=((yreal.l & 0x7f000000) >> 24)-64;
 mantisse=(yreal.l & 0x00ffffff);
 if (vorzeichen)
  return(-mantisse/16777216.0*pow(16.0,(double)exponent));
 else
  return( mantisse/16777216.0*pow(16.0,(double)exponent));
}


/* read one complete SEGY "tape block" (with block labels) */
long SEGYread (char *buffer, int  size, FILE *fp)
{
 long block1, block2;

 if (fread((char *)&block1, 1, 4, fp)==4)
 {
  if (block1 > size)
  {
   fprintf(stderr, "Block too large ! \n");
   fprintf(stderr, "given : %d(%08x) < found : %d(%08x) --- exit\n",
                   size, size, block1, block1);
   exit(98);
  }
  fread(buffer, 1, block1, fp);
  fread((char *)&block2, 1, 4, fp);
  if (block2 != block1)
  {
   fprintf(stderr, "Wrong block structure ! \n");
   fprintf(stderr, "start : %d(%08x)  != end : %d(%08x) --- exit\n",
                   block1, block1, block2, block2);
   exit(99);
  }
  return(block1);
 }
 else /* eof */
  return(0);
}

/* write one complete SEGY "tape block" (with block labels) */
void SEGYwrite(char *buffer, int  size, FILE *fp)
{
 long block;

 block=size;
 fwrite((char *)&block, 1, 4, fp);
 fwrite(buffer, 1, size, fp);
 fwrite((char *)&block, 1, 4, fp);
}

FILE *segy_open_read(char *name,struct asciiheader *ascii,struct binaryheader *binary)
{
 FILE *fp;

 fp=fopen(name, "r");
 if (fp!=NULL)
 {
  if (SEGYread((char *)ascii, 3200, fp)!=3200)
/* unexpected eof */
  {
   fclose(fp);
   return(NULL);
  }
  if (SEGYread((char *)binary, 400, fp)!=400)
/* unexpected eof */
  {
   fclose(fp);
   return(NULL);
  }
 }
 return(fp);
}

FILE *segy_open_write(char *name,struct asciiheader *ascii,struct binaryheader *binary)
{
 FILE *fp;
 struct binaryheader bincopy;

 fp=fopen(name, "w");
 if (fp!=NULL)
 {
  SEGYwrite((char *)ascii, 3200, fp);
  memcpy(&bincopy, binary, sizeof(struct binaryheader));
  bincopy.format=1;
  SEGYwrite((char *)&bincopy, 400, fp);
 }
 return(fp);
}

long segy_read(FILE *fp,struct traceheader *trace,float *data,struct binaryheader *binary)
{
 long i;

 i=SEGYread((char *)&work,
            (int)(240+(binary->format==3 ? 2 : 4)*binary->lproc), fp);
 if (i<240)
/* incomplete trace (eof) */
 {
  trace->lproc=0;
  return(0);
 }
 memcpy((char *)trace, (char *)&work.th, 240);
 if (i!=240+(binary->format==3 ? 2 : 4)*trace->lproc)
/* wrong sample count */
 {
  trace->lproc=0;
  return(0);
 }
 for (i=0; i<trace->lproc; i++)
/* Daten umrechnen */
  switch (binary->format)
  {
   case 1 : data[i]=real2float((char *)&work.da.l[i]);
                                               break;
   case 2 : data[i]=(float)work.da.l[i];       break;
   case 3 : data[i]=(float)work.da.s[i];       break;
   case 4 : data[i]=(float)(work.da.g[i].v*(1<<work.da.g[i].e));
                                               break;
  }
 return((long)trace->lproc);
}

long segy_write(FILE *fp,struct traceheader *trace, float *data)
{
 long i;

 memcpy((char *)&work.th, (char *)trace, 240);
 for (i=0; i<trace->lproc; i++)
  float2real(data[i],(char *)&work.da.l[i]);
 SEGYwrite((char *)&work, (int)(240+4*trace->lproc), fp);
 return((long)trace->lproc);
}

FILE *image_open_write(char *name,struct asciiheader *ascii,struct binaryheader *binary)
{
 FILE *fp;
 struct binaryheader bincopy;

 fp=fopen(name, "w");
 if (fp!=NULL)
 {
  write((char *)ascii, 1, 3200, fp);
  memcpy(&bincopy, binary, sizeof(struct binaryheader));
  bincopy.format=1;
  write((char *)&bincopy, 1, 400, fp);
 }
 return(fp);
}

long image_write(FILE *fp,struct traceheader *trace,float *data)
{
 long i;

 memcpy((char *)&work.th, (char *)trace, 240);
 for (i=0; i<trace->lproc; i++)
  float2real(data[i],(char *)&work.da.l[i]);
 write((char *)&work, 1, (int)(240+4*trace->lproc), fp);
 return((long)trace->lproc);
}