segdread.c

su 的源代码库

  */

static void F8038_to_float (Sfio_t *from, float to[], int len)
{
 int i;
 long int ex1_4;
 long int ex2_4;
 long int value;

 for (i = 0; i < len; i ++) {
      ex1_4 = GET_S(from);
      ex2_4 = GET_S(from);
      value = (ex1_4<<16) | (ex2_4&65535);
      *(to++) = (float) value;
 }
} 

/* F8042_to_float - convert 8 bit hexadecimal demultiplexed data into floating numbers
 *
 * Credits:
 *      SEP:  Stew Levin
 *
 * Parameters:
 *    from   - input sfio unit
 *    to     - output vector
 *    len    - number of floats in vector
 *
 */
/*
 *
 * Format 8042 is a 1 byte per word representation.
 * According to the SEG specifications, the bit
 * layout of the byte is:
 *
 *
 *  Bit       0     1     2     3     4     5     6     7
 *-----------------------------------------------------------
 * Byte 1     S    C1    C0    Q-1   Q-2   Q-3   Q-4   Q-5
 *
 * S=sign bit. - (One = negative number)
 * C=hexadecimal exponent. - This is a 2 bit positive binary exponent of 16
 *                CC
 *   written as 16    where CC can assume values of 0-3.
 * Q1-5-fraction. - This is a 5 bit positive binary fraction.
 *   The radix point is to the left of the most significant bit (Q-1)
 *                                  -1
 *   with the MSB being defined as 2 .  The fraction can have values
 *           -5        -5
 *   from 1-2   to -1+2  .  
 *                             CC    MP                   MP
 * Input signal = S.QQQQ,Q x 16   x 2   millivolts where 2    is the
 *   value required to descale the data word to the recording system
 *   input level.  MP is defined in Byte 8 of each of the corre-
 *   sponding channel set descriptors in the scan type header.
 */

static void F8042_to_float (Sfio_t *from, float to[], int len)
{
 register int i;
 register int ex1_4;
 int expo;
 short fraction;

 for (i = 0; i < len; i ++) {
      ex1_4 = GET_C(from);
      expo = ((ex1_4 >> 3) & 12) - 5;
      fraction = ex1_4 & 31;
      if (ex1_4 & 128) fraction = -fraction;
      *(to++) = ldexp((double) fraction, expo);
 }
} 

/* F8044_to_float - convert 16 bit hexadecimal demultiplexed data into floating numbers
 *
 * Credits:
 *      SEP:  Stew Levin
 *
 * Parameters:
 *    from   - input sfio unit
 *    to     - output vector
 *    len    - number of floats in vector
 *
 */
/*
 *
 * Format 8044 is a 2 byte per word representation.
 * According to the SEG specifications, the bit
 * layout of the bytes is:
 *
 *
 *  Bit       0     1     2     3     4     5     6     7
 *-----------------------------------------------------------
 * Byte 1     S    C1    C0    Q-1   Q-2   Q-3   Q-4   Q-5
 * Byte 2    Q-6   Q-7   Q-8   Q-9   Q-10  Q-11  Q-12  Q-13
 *
 * S=sign bit. - (One = negative number)
 * C=hexadecimal exponent. - This is a 2 bit positive binary exponent of 16
 *                CC
 *   written as 16    where CC can assume values of 0-3.
 * Q1-13-fraction. - This is a 13 bit positive binary fraction.
 *   The radix point is to the left of the most significant bit (Q-1)
 *                                  -1
 *   with the MSB being defined as 2 .  The fraction can have values
 *           -13        -13
 *   from 1-2    to -1+2   .  
 *                                       CC    MP                   MP
 * Input signal = S.QQQQ,QQQQ,QQQQ,Q x 16   x 2   millivolts where 2    
 *   is the value required to descale the data word to the recording system
 *   input level.  MP is defined in Byte 8 of each of the corre-
 *   sponding channel set descriptors in the scan type header.
 */

static void F8044_to_float (Sfio_t *from, float to[], int len)
{
 register int i;
 register int ex1_4;
 int expo;
 short fraction;

 for (i = 0; i < len; i ++) {
      ex1_4 = GET_S(from);
      expo = ((ex1_4 >> 11) & 12) - 13;
      fraction = ex1_4 & 8191;
      if (ex1_4 & 32768) fraction = -fraction;
      *(to++) = ldexp((double) fraction, expo);
 }
} 

/* F8048_to_float - convert 32 bit hexadecimal demultiplexed data into floating numbers
 *
 * Credits:
 *      SEP:  Stew Levin
 *
 * Parameters:
 *    from   - input sfio unit
 *    to     - output vector
 *    len    - number of floats in vector
 *
 */
/*
 *
 * Format 8048 is a 4 byte per word representation.
 * According to the SEG specifications, the bit
 * layout of the bytes is:
 *
 *
 *  Bit       0     1     2     3     4     5     6     7
 *-----------------------------------------------------------
 * Byte 1     S    C6    C5    C4    C3    C2    C1    C0 
 * Byte 2    Q-1   Q-2   Q-3   Q-4   Q-5   Q-6   Q-7   Q-8 
 * Byte 3    Q-9   Q-10  Q-11  Q-12  Q-13  Q-14  Q-15  Q-16
 * Byte 4    Q-17  Q-18  Q-19  Q-20  Q-21  Q-22  Q-23  0
 *
 * S=sign bit. - (One = negative number)
 * C=hexadecimal exponent. - This is a binary exponent of 16
 *                (CCCCCCC-64)
 *   written as 16             where CC can assume values of 0-127.
 * Q1-23-fraction. - This is a 23 bit positive binary fraction.
 *   The radix point is to the left of the most significant bit (Q-1)
 *                                  -1
 *   with the MSB being defined as 2 .  The sign and fraction can have
 *                 -23        -23
 *   values from 1-2    to -1+2   .  
 *                                   C-64    MP                   MP
 * Input signal = S.QQQQ,...,QQQ x 16     x 2   millivolts where 2    
 *   is the value required to descale the data word to the recording system
 *   input level.  MP is defined in Byte 8 of each of the corre-
 *   sponding channel set descriptors in the scan type header.
 *   The data recording method has more than sufficient range to
 *   handle the dynamic range of a typical seismic system.  Thus, MP
 *   may not be needed to account for any scaling and may be recorded
 *   as zero.
 */

static void F8048_to_float (Sfio_t *from, float to[], int len)
{
 register int i;
 register int ex1_4;
 int expo;
 long int fraction;

 for (i = 0; i < len; i ++) {
      ex1_4 = GET_S(from); 
      expo = ((ex1_4 >> 6) & 508) - (24+256);
      fraction = ex1_4 & 255;
      fraction <<= 16; fraction |= (GET_S(from)&65535);
      if (ex1_4 & 32768) fraction = -fraction;
      *(to++) = ldexp((double) fraction, expo);
 }
} 

/* F8058_to_float - convert 32 bit IEEE float demultiplexed data into floating numbers
 *
 * Credits:
 *      SEP:  Stew Levin
 *
 * Parameters:
 *    from   - input sfio unit
 *    to     - output vector
 *    len    - number of floats in vector
 *
 */
/*
 *
 * Format 8058 is a 4 byte per word representation.
 * According to the SEG specifications, the bit
 * layout of the bytes is:
 *
 *
 *  Bit       0     1     2     3     4     5     6     7
 *-----------------------------------------------------------
 * Byte 1     S    C7    C6    C5    C4    C3    C2    C1 
 * Byte 2    C0    Q-1   Q-2   Q-3   Q-4   Q-5   Q-6   Q-7
 * Byte 3    Q-8   Q-9   Q-10  Q-11  Q-12  Q-13  Q-14  Q-15
 * Byte 4    Q-16  Q-17  Q-18  Q-19  Q-20  Q-21  Q-22  Q-23
 *
 * S=sign bit. - (One = negative number)
 * C=exponent. - This is a excess-127 binary exponent of 2
 *               (CCCCCCCC-127)
 *   written as 2               where CC can assume values of 0-255.
 * Q1-23-fraction. - This is a 23 bit positive binary fraction.
 *   The radix point is to the left of the most significant bit (Q-1)
 *                                  -1
 *   with the MSB being defined as 2 .  With the exceptions noted below:
 *
 *                    S                    C-127    MP                   MP
 * Input signal = (-1) x 1.QQQQ,...,QQQ x 2     x 2   millivolts where 2    
 *   is the value required to descale the data word to the recording system
 *   input level.  MP is defined in Byte 8 of each of the corre-
 *   sponding channel set descriptors in the scan type header.
 *   The data recording method has more than sufficient range to
 *   handle the dynamic range of a typical seismic system.  Thus, MP
 *   may not be needed to account for any scaling and may be recorded
 *   as zero.
 *
 * Exceptions:
 *
 * If C=0 then
 *                    S                    -126     MP
 * Input signal = (-1) x 0.QQQQ,...,QQQ x 2     x 2   millivolts
 *
 * If C=255 and Q=0, then
 *                    S 
 * Input signal = (-1) x infinity  (overflow)
 *
 * If C=255 and Q!=0, then
 *                      
 * Input signal = NaN  (Not-a-Number)
 */

static void F8058_to_float (Sfio_t *from, float to[], int len)
{
 register int i;
 register int ex1_4, ex2_4;
 int expo;
 long int fraction;

 for (i = 0; i < len; i ++) {
      ex1_4 = GET_S(from); 
      ex2_4 = GET_S(from);
      expo = ((ex1_4 >> 7) & 255);
      fraction = ex1_4 & 127;
      fraction <<= 16; fraction |= (ex2_4&65535);
      if(expo) fraction |= 8388608;
      else fraction <<= 1;
      if (ex1_4 & 32768) fraction = -fraction;
      *(to++) = ldexp((double) fraction, expo-(23+127));
 }
} 

int
get_gh1(general_header_1 * gh1, Sfio_t *tapeun)
{
  int status;

  gh1->f[0]=GET_UC(tapeun);
  gh1->f[1]=GET_UC(tapeun);
  gh1->y=GET_US(tapeun);
  gh1->k1_k2=GET_C(tapeun);
  gh1->k3_k4=GET_C(tapeun);
  gh1->k5_k6=GET_C(tapeun);
  gh1->k7_k8=GET_C(tapeun);
  gh1->k9_k10=GET_C(tapeun);
  gh1->k11_k12=GET_C(tapeun);
  gh1->yr=GET_UC(tapeun);
  gh1->gh_dy1=GET_UC(tapeun);
  
  gh1->dy=GET_UC(tapeun);
  gh1->h=GET_UC(tapeun);
  gh1->mi=GET_UC(tapeun);
  gh1->se=GET_UC(tapeun);
  
  gh1->m[0]=GET_UC(tapeun);
  gh1->m[1]=GET_UC(tapeun);
  gh1->m[2]=GET_UC(tapeun);
  gh1->b[0]=GET_UC(tapeun);
  
  gh1->b[1]=GET_UC(tapeun);
  gh1->b[2]=GET_UC(tapeun);
  gh1->i=GET_UC(tapeun);
  gh1->p_sbx=GET_UC(tapeun);
  
  gh1->sb=GET_UC(tapeun);
  gh1->z_r1=GET_UC(tapeun);
  gh1->r=GET_UC(tapeun);
  gh1->str=GET_UC(tapeun);
  
  gh1->cs=GET_UC(tapeun);
  gh1->sk=GET_UC(tapeun);
  gh1->ec=GET_UC(tapeun);
  gh1->ex=GET_UC(tapeun);

  status = (sferror(tapeun)||sfeof(tapeun));
  
  return (status?EXIT_FAILURE:EXIT_SUCCESS);
}

int
get_gh2(general_header_2 * gh2, Sfio_t *tapeun)
{
  int status;
  
  gh2->ef[0]=GET_UC(tapeun);
  gh2->ef[1]=GET_UC(tapeun);
  gh2->ef[2]=GET_UC(tapeun);
  gh2->en[0]=GET_UC(tapeun);
  
  gh2->en[1]=GET_UC(tapeun);
  gh2->ecx[0]=GET_UC(tapeun);
  gh2->ecx[1]=GET_UC(tapeun);
  gh2->eh[0]=GET_UC(tapeun);
  
  gh2->eh[1]=GET_UC(tapeun);
  gh2->x1=GET_C(tapeun);
  gh2->rev[0]=GET_UC(tapeun);
  gh2->rev[1]=GET_UC(tapeun);
  
  gh2->gt=GET_US(tapeun);
  gh2->erl[0]=GET_UC(tapeun);
  gh2->erl[1]=GET_UC(tapeun);
  
  gh2->erl[2]=GET_UC(tapeun);
  gh2->x2=GET_C(tapeun);
  gh2->bn=GET_UC(tapeun);
  gh2->x3[0]=GET_C(tapeun);
  
  gh2->x3[1]=GET_C(tapeun);
  gh2->x3[2]=GET_C(tapeun);
  gh2->x3[3]=GET_C(tapeun);
  gh2->x3[4]=GET_C(tapeun);
  
  gh2->x3[5]=GET_C(tapeun);
  gh2->x3[6]=GET_C(tapeun);
  gh2->x3[7]=GET_C(tapeun);
  gh2->x3[8]=GET_C(tapeun);
  
  gh2->x3[9]=GET_C(tapeun);
  gh2->x3[10]=GET_C(tapeun);
  gh2->x3[11]=GET_C(tapeun);
  gh2->x3[12]=GET_C(tapeun);

  status = (sferror(tapeun)||sfeof(tapeun));

  return (status?EXIT_FAILURE:EXIT_SUCCESS);
}

int
get_ghn(general_header_n * ghn, Sfio_t *tapeun)
{
  int status;
  
  ghn->x1[0]=GET_C(tapeun);
  ghn->x1[1]=GET_C(tapeun);