Using SAC

Using the SAC Library

Overview

In addition to being able to read and write SAC data files in one's own C or FORTRAN programs (see SAC Reading and Writing Routines ), one can use many of SAC's data-processing routines in stand-alone codes if one uses specific flags in the compiling stage and the SAC libraries in the linking stage. These libraries libsac.a and libsacio.a are in ${SACHOME}/lib.

If ${SACHOME} is /usr/local/sac, the following commands will compile the Fortran anc C programs for filter.f and filterc.c respectively:

f77 -o filterf filter.f -I/usr/local/sac/include  -L/usr/local/sac/lib  -lsacio -lsac

gcc -o filterc filterc.c -I/usr/local/sac/include  -L/usr/local/sac/lib  -lsacio -lsac

Filtering - Fortran

      program filter
      implicit none

      include "sacf.h"

!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=1000)

!     Define the Data Array of size MAX
      real yarray, xarray
      dimension yarray(MAX)

!     Declare Variables used in the rsac1() subroutine
      real beg, delta
      integer nlen
      character*20 KNAME
      integer nerr

!     Define variables used in the filtering routine
      real *8 low, high
      real *8 transition_bandwidth, attenuation
      real *8 delta_d
      integer order, passes

      kname = 'filterf_in.sac'

      call rsac1(kname, yarray, nlen, beg, delta, MAX, nerr)
      delta_d = delta

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

      low    = 0.10
      high   = 1.00
      passes = 2
      order  = 4
      transition_bandwidth = 0.0
      attenuation = 0.0
!     Call xapiir ( Apply a IIR Filter )
!        - yarray - Original Data
!        - nlen   - Number of points in yarray
!        - proto  - Prototype of Filter
!                 - SAC_FILTER_BUTTERWORK        - Butterworth
!                 - SAC_FILTER_BESSEL            - Bessel
!                 - SAC_FILTER_CHEBYSHEV_TYPE_I  - Chebyshev Type I
!                 - SAC_FILTER_CHEBYSHEV_TYPE_II - Chebyshev Type II
!        - transition_bandwidth (Only for Chebyshev Filter)
!                 - Bandwidth as a fraction of the lowpass prototype
!                   cutoff frequency
!        - attenuation (Only for Chebyshev Filter)
!                 - Attenuation factor, equals amplitude reached at
!                   stopband egde
!        - order  - Number of poles or order of the analog prototype
!                   4 - 5 should be ample
!                   Cannot exceed 10
!        - type   - Type of Filter
!                 - SAC_FILTER_BANDPASS
!                 - SAC_FILTER_BANDREJECT
!                 - SAC_FILTER_LOWPASS
!                 - SAC_FILTER_HIGHPASS
!        - low    - Low Frequency Cutoff [ Hertz ]
!                   Ignored on SAC_FILTER_LOWPASS
!        - high   - High Frequency Cutoff [ Hertz ]
!                   Ignored on SAC_FILTER_HIGHPASS
!        - delta  - Sampling Interval [ seconds ]
!        - passes - Number of passes 
!                 - 1 Forward filter only
!                 - 2 Forward and reverse (i.e. zero-phase) filtering
!
      call xapiir(yarray, nlen, 
     +            SAC_BUTTERWORTH, 
     +            transition_bandwidth, attenuation,
     +            order, 
     +            SAC_BANDPASS,
     +            low, high, delta_d, passes)

!     Do more processing ....
      
      xarray = 0
      kname='filterf_out.sac'
!     Write the SAC file kname
!       - kname holds the name of the file to be written
!       - xdata Input Time Data      
!       - yfunc Input Amplitude Data
!       - nerr Error return Flag
      call wsac0(kname,xarray,yarray,nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program filter

Filtering - C


#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "sac.h"
#include "sacio.h"

#define MAX 1000

int 
main(int argc, char *argv[]) {

    /* Local variables */
    double low, high, attenuation, transition_bandwidth;;
    int nlen, nerr, max;

    float beg, delta;
    double delta_d;
    char *kname;
    int order;

    int passes;
    float yarray[1000];
    float xarray[1];

    max = MAX;

    /*     Define the Maximum size of the data Array 
     * Filter Prototypes 
     * Filter Types 
     *     Define the Data Array of size MAX 
     *     Declare Variables used in the rsac1() subroutine 
     *     Define variables used in the filtering routine 
     */
        
    kname = strdup("filterc_in.sac");
    rsac1(kname, yarray, &nlen, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);
    delta_d = delta;
    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    low    = 0.10;
    high   = 1.00;
    passes = 2;
    order  = 4;
    transition_bandwidth = 0.0;
    attenuation = 0.0;

    /*     Call xapiir ( Apply a IIR Filter ) 
     *        - yarray - Original Data 
     *        - nlen   - Number of points in yarray 
     *        - proto  - Prototype of Filter 
     *                 - SAC_FILTER_BUTTERWORK        - Butterworth 
     *                 - SAC_FILTER_BESSEL            - Bessel 
     *                 - SAC_FILTER_CHEBYSHEV_TYPE_I  - Chebyshev Type I 
     *                 - SAC_FILTER_CHEBYSHEV_TYPE_II - Chebyshev Type II 
     *        - transition_bandwidth (Only for Chebyshev Filter) 
     *                 - Bandwidth as a fraction of the lowpass prototype 
     *                   cutoff frequency 
     *        - attenuation (Only for Chebyshev Filter) 
     *                 - Attenuation factor, equals amplitude reached at 
     *                   stopband egde 
     *        - order  - Number of poles or order of the analog prototype 
     *                   4 - 5 should be ample 
     *                   Cannot exceed 10 
     *        - type   - Type of Filter 
     *                 - SAC_FILTER_BANDPASS 
     *                 - SAC_FILTER_BANDREJECT 
     *                 - SAC_FILTER_LOWPASS 
     *                 - SAC_FILTER_HIGHPASS 
     *        - low    - Low Frequency Cutoff [ Hertz ] 
     *                   Ignored on SAC_FILTER_LOWPASS 
     *        - high   - High Frequency Cutoff [ Hertz ] 
     *                   Ignored on SAC_FILTER_HIGHPASS 
     *        - delta  - Sampling Interval [ seconds ] 
     *        - passes - Number of passes 
     *                 - 1 Forward filter only 
     *                 - 2 Forward and reverse (i.e. zero-phase) filtering 
     */
    xapiir(yarray, nlen, SAC_BUTTERWORTH, 
           transition_bandwidth, attenuation, 
           order, 
           SAC_BANDPASS, 
           low, high, 
           delta_d, passes);

    /*     Do more processing ....  */
    xarray[0] = 0;
    kname = strdup("filterc_out.sac");    
    wsac0(kname, xarray, yarray, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }
    

    return 0;
}

Envelope Function - Fortran

      program envelopef
      implicit none

      include "sacf.h"

!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=1000)

!     Define the Data Array of size MAX
      real yarray, xarray, yenv
      dimension yarray(MAX), yenv(MAX)

!     Declare Variables used in the rsac1() subroutine
      real beg, delta
      integer nlen
      character*20 KNAME
      integer nerr

      kname = 'envelopef_in.sac'

      call rsac1(kname, yarray, nlen, beg, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Call envelope ( Envelope Routine )
!        - nlen   - Number of points in yarray
!        - yarray - Original Data
!        - yenv   - Envelope of the Original Data
!
      call envelope(nlen, yarray, yenv)

!     Do more processing ....
      
      xarray = 0
      kname='envelopef_out.sac'
!     Write the SAC file kname
!       - kname holds the name of the file to be written
!       - xdata Input Time Data      
!       - yfunc Input Amplitude Data
!       - nerr Error return Flag
      call wsac0(kname,xarray,yenv,nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program envelopef

Envelope Function - C


#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include "sac.h"
#include "sacio.h"

#define MAX 1000

int 
main(int argc, char *argv[]) {

    /* Local variables */
    int nlen, nerr, max;

    float beg, delta;
    char *kname;

    float yarray[1000];
    float xarray[1];
    float *yenv;

    max = MAX;

    /*     Define the Maximum size of the data Array 
     *     Define the Data Array of size MAX 
     *     Declare Variables used in the rsac1() subroutine 
     *     Define variables used in the envelope routine 
     */
        
    kname = strdup("envelopec_in.sac");
    rsac1(kname, yarray, &nlen, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    /* Allocate space for the envelope of yarray */
    yenv = (float *) malloc(sizeof(float) * nlen);
    if(yenv == NULL) {
      fprintf(stderr, "Error allocating memory for envelope\n");
      exit(-1);
    }

    /*     Call envelope ( Envelope Function )
     *        - nlen   - Number of points in yarray
     *        - yarray - Input array to take the envelope of
     *        - yenv   - Envelope of yarray
     */
    envelope(nlen, yarray, yenv);

    /*     Do more processing ....  */
    xarray[0] = 0;
    kname = strdup("envelopec_out.sac");    
    wsac0(kname, xarray, yenv, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }
    
    free(yenv);

    return 0;
}

Correlation - Fortran

      program envelopef
      implicit none

      include "sacf.h"

      integer i
!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=4000)

!     Define the Data Array of size MAX
      real yarray1, yarray2, xarray, out, ytmp
      dimension yarray1(MAX), yarray2(MAX), out(MAX*4), ytmp(MAX*4)

!     Declare Variables used in the rsac1() subroutine
      real beg1, beg2, delta, endv
      integer nlen, nlen1, nlen2
      character*30 KNAME
      integer nerr
      real fval
      
!     Cross Correlation Variables
      integer nwin, wlen, nfft
      character *256 error
      real max_value, max_time

!     Read in the first data file
      kname = 'correlatef_in1.sac'
      call rsac1(kname, yarray1, nlen1, beg1, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Read in the second data file
      kname = 'correlatef_in2.sac'
      call rsac1(kname, yarray2, nlen2, beg2, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

      nlen = nlen1
!
!     If the signals are not the same length, then find the longest
!     signal, make both signals that length by filling the remainder
!     with zeros (pad at the end) and then run them through crscor
!     This should be fixed in upcoming releases and the introduction
!     of a "correlate" function so you do not need to handle the 
!     signal length, padding, and window lengths, ...
!

      nwin = 1
      wlen = nlen
      nfft = 0
!     Call crscor ( Cross Correlation )
!        - yarray1 - First  Input array to correlate
!        - yarray2 - Second Input array to correlate
!        - nlen    - Number of points in yarray and yarray2
!        - nwin    - Windows to use in the correlation
!        - wlen    - Length of the windows
!        - type    - Type of Window (SAC_RECTANGLE)
!        - out     - output sequence 
!        - nfft    - Length of the output sequence
!        - error   - Error Message
!

      call crscor(yarray1, yarray2, nlen, 
     &            nwin, wlen, SAC_RECTANGLE,
     &            ytmp, nfft, error)

      do i = 1,MAX*4
         out(i) = 0.0
      enddo
!     
!     out[1 : nlen1 - 1 ] <-- ytmp[ nfft - nlen1 + 2 : nfft  ]
!     out[nlen1 : nlen1 + nlen2 - 1 ] <-- ytmp[ 1 : nlen2 ]
!
      do i = 1,nlen1-1
         out(i) = ytmp(nfft - nlen1 + i + 1)
      enddo
      do i = 1,nlen2
         out(nlen1 + i - 1) = ytmp(i)
      enddo


      nfft = nlen1 + nlen2 - 1
      xarray = 0
      beg1 = -delta * (nlen1 - 1.0) + (beg2 - beg1)
      endv = beg1 + delta * (nfft - 1)

      call setnhv('npts',   nfft,    nerr)
      call setfhv('delta',  delta,   nerr)
      call setlhv('leven',  .true.,  nerr)
      call setfhv('b',      beg1,    nerr)
      call setfhv('e',      endv,    nerr)
      call setihv('iftype', "itime", nerr)
      call setnhv('nzyear', SAC_NUMBER_UNDEFINED, nerr)
      call setnhv('nzhour', SAC_NUMBER_UNDEFINED, nerr)
!     Find the maximum value and time of the correlation function 
      max_value = out(1)
      max_time  = 0
      do i = 2,nfft
         if(out(i) > max_value) then
            max_value = out(i)
!           Negative shifts are at the end of the correlation sequence
            if(i > nfft/2) then
               max_time = beg1 + (i - nfft) * delta
            else
               max_time = beg1 + i * delta
            endif
         endif
      enddo

!      call setfhv( 'user0', max_time,  nerr)
!      call setfhv( 'user1', max_value, nerr)


!     Write the SAC file
      kname='correlatef_out1.sac'
      call wsac0(kname, xarray, out, nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname
          call exit(-1)
      endif

      call exit(0)

      end program envelopef

Correlation - C


#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <sac.h>
#include <sacio.h>

#define MAX        4000
#define ERROR_MAX  256

int 
main(int argc, char *argv[]) {

    /* Local variables */
    int i;
    int nlen, nlen1, nlen2, nerr, max;

    float beg1, beg2, delta, end;
    char *kname;

    float yarray1[MAX], yarray2[MAX], ytmp[MAX*4], xarray[1];
    float *out;
    float max_value, max_time;

    int nwin, wlen, nfft, leven, when;

    char error[ERROR_MAX];

    max = MAX;

    /* Read in the first file  */        
    kname = strdup("correlatec_in1.sac");
    rsac1(kname, yarray1, &nlen1, &beg1, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file(%d): %s\n", nerr, kname);
      exit(-1);
    }

    /* Read in the second file  */
    kname = strdup("correlatec_in2.sac");
    rsac1(kname, yarray2, &nlen2, &beg2, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    nlen = nlen1;
    /** 
     *  If the signals are not the same length, then find the longest
     *  signal, make both signals that length by filling the remainder
     *  with zeros (pad at the end) and then run them through crscor
     *  This should be fixed in upcoming releases and the introduction
     *  of a "correlate" function so you do not need to handle the 
     *  signal length, padding, and window lengths, ...
     *
     */

    /* Allocate space for the correlation of yarray1 and yarray2 */
    max = next2((2 * nlen) - 1) * 2;
    out = (float *) malloc(sizeof(float) * max);
    if(out == NULL) {
      fprintf(stderr, "Error allocating memory for correlation\n");
      exit(-1);
    }

    /* Set up values for the cross correlation */
    nwin = 1;
    wlen = nlen;
    nfft = 0;
    
    /*     Call crscor ( Cross Correlation )
     *        - yarray1 - First  Input array to correlate
     *        - yarray2 - Second Input array to correlate
     *        - nlen    - Number of points in yarray and yarray2
     *        - nwin    - Windows to use in the correlation
     *        - wlen    - Length of the windows
     *        - type    - Type of Window (SAC_RECTANGLE)
     *        - out     - output sequence 
     *        - nfft    - Length of the output sequence
     *        - error   - Error Message
     *        - err_len - Length of Error Message (on input)
     */
    crscor(yarray1, yarray2, nlen, 
           nwin, wlen, SAC_RECTANGLE,
           ytmp, &nfft, error, ERROR_MAX);
    
    for(i = 0; i < max; i++) {
      out[i] = 0;
    }
    
    /*     
     *     out[0 : nlen1 - 2 ] <-- ytmp[ nfft - nlen1 + 1 : nfft -1 ]
     *     out[nlen1 - 1 : nlen1 + nlen2 - 2 ] <-- ytmp[ 0 : nlen2-1 ]
     */
    for(i = 0; i <= nlen1 - 2; i++) {
      out[i] = ytmp[nfft - nlen1 + i + 1];
    }
    for(i = 0; i <= nlen2 - 1; i++) {
      out[nlen1 + i - 1] = ytmp[i];
    }


    nfft = nlen1 + nlen2 - 1;
    xarray[0] = 0;
    leven = TRUE;
    beg1 = -delta * (nlen1 - 1) + (beg2 - beg1);
    end = beg1 + delta * (nfft - 1);

    setnhv ( "npts",   &nfft,   &nerr, SAC_STRING_LENGTH);
    setfhv ( "delta",  &delta,  &nerr, SAC_STRING_LENGTH);
    setlhv ( "leven",  &leven,  &nerr, SAC_STRING_LENGTH);
    setfhv ( "b",      &beg1,   &nerr, SAC_STRING_LENGTH);
    setfhv ( "e",      &end,    &nerr, SAC_STRING_LENGTH);
    setihv ( "iftype", "itime", &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);
    when = SAC_NUMBER_UNDEFINED;
    setnhv ( "nzyear", &when,   &nerr, SAC_STRING_LENGTH);
    setnhv ( "nzhour", &when,   &nerr, SAC_STRING_LENGTH);

    /* Find the maximum value and time of the correlation function */
    max_value = out[0];
    max_time  = 0;
    for(i = 1; i < nfft; i++) {
      if(out[i] > max_value) {
        max_value = out[i];
        /* Negative shifts are at the end of the correlation sequence */
        if(i > nfft/2) { 
          max_time  = (i - nfft) * delta;
        } else {
          max_time  = i * delta;
        }
      }
    }

    /*
      setfhv( "user0", &max_time,  &nerr, SAC_STRING_LENGTH);
      setfhv( "user1", &max_value, &nerr, SAC_STRING_LENGTH);
    */

    /*   Write out the correlation function   */
    kname = strdup("correlatec_out1.sac");    
    wsac0(kname, xarray, out, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }
    
    free(out);

    return 0;
}

Convolution - Fortran

      program envelopef
      implicit none

      include "sacf.h"

      integer i,j
!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=4000)

!     Define the Data Array of size MAX
      real yarray1, yarray2, ytmp, xarray, out
      dimension yarray1(MAX), yarray2(MAX), ytmp(MAX), out(MAX*4)

!     Declare Variables used in the rsac1() subroutine
      real beg, delta, endv
      integer nlen1, nlen2, nlen
      character*30 KNAME
      integer nerr
      
!     Cross Correlation Variables
      integer nwin, wlen, nfft
      character *256 error
      character *24  kevnm

!     Read in the first data file
      kname = 'convolvef_in1.sac'
      call rsac1(kname, ytmp, nlen1, beg, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Read in the second data file
      kname = 'convolvef_in2.sac'
      call rsac1(kname, yarray2, nlen2, beg, delta, MAX, nerr)

      if(nerr .NE. 0) then
          write(*,*)'Error reading in file: ',kname
          call exit(-1)
      endif

!     Reverse the First Signal */      
      j = 1
      do i = nlen1,1,-1
         yarray1(j) = ytmp(i)
         j = j + 1
      enddo

      nlen = nlen1
      if(nlen2 > nlen) then
         nlen = nlen2
      endif


      nwin = 1
      wlen = nlen
      nfft = 0
!     Call crscor ( Cross Correlation )
!        - yarray1 - First  Input array to correlate
!        - yarray2 - Second Input array to correlate
!        - nlen    - Number of points in yarray and yarray2
!        - nwin    - Windows to use in the correlation
!        - wlen    - Length of the windows
!        - type    - Type of Window (SAC_RECTANGLE)
!        - out     - output sequence 
!        - nfft    - Length of the output sequence
!        - error   - Error Message
!

      call crscor(yarray1, yarray2, nlen, 
     &            nwin, wlen, SAC_RECTANGLE,
     &            out, nfft, error)


!     Zero out the tmp signal      
      do i = 1, MAX
         ytmp(i) = 0.0
      enddo

!     Reconstruct the signal from the "cross correlation" back to front
!  
!     ytmp[1        : nlen1 - 1         ] <- out[nfft-nlen1+2 : nfft  ] 
!     ytmp[nlen1    : nlen1 + nlen2 - 1 ] <- out[1            : nlen2 ]
!
!     nfft is the last point of the output sequence
!
      do i = 1,nlen1-1
         ytmp(i) = out(nfft - nlen1 + i + 1)
      enddo
      do i = 1, nlen2
         ytmp(nlen1 + i - 1) = out(i)
      enddo


      nfft = nlen1 + nlen2 - 1
      xarray = 0
      beg = 0
      endv = beg + delta * (nfft - 1)
      j = 1

      call newhdr()
      call setnhv('npts',   nfft,    nerr)
      call setfhv('delta',  delta,   nerr)
      call setlhv('leven',  .true.,  nerr)
      call setfhv('b',      beg,     nerr)
      call setfhv('e',      endv,    nerr)
      call setihv('iftype', 'itime', nerr)
      call setkhv('kstnm',  'sta',   nerr)
      call setkhv('kcmpnm', 'Q',     nerr)
      call setnhv('nwfid',  j, nerr)
      kevnm = 'FUNCGEN: TRIANGLE'
      call setkhv ('kevnm', kevnm, nerr)
!     Write the SAC file
      kname='convolvef_out1.sac'
      call wsac0(kname, xarray, ytmp, nerr)
      if(nerr .NE. 0) then
          write(*,*)'Error writing out file: ',kname,nerr
          call exit(-1)
      endif

      call exit(0)

      end program envelopef

Convolution - C


#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <sac.h>
#include <sacio.h>

#define MAX        4000
#define ERROR_MAX  256

int 
main(int argc, char *argv[]) {

    /* Local variables */
    int i, j;
    int nlen, nlen1, nlen2, nerr, max;

    float beg, delta, end;
    char *kname;

    float yarray1[MAX], yarray2[MAX], ytmp[MAX], xarray[1];
    float *out;

    int nwin, wlen, nfft, leven;

    char error[ERROR_MAX];

    max = MAX;

    for(i = 0; i < MAX; i++) {
      yarray1[i] = 0.0;
      yarray2[i] = 0.0;
      ytmp[i] = 0.0;
    }
    /* Read in the first file  */        
    kname = strdup("convolvec_in1.sac");
    rsac1(kname, ytmp, &nlen1, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file(%d): %s\n", nerr, kname);
      exit(-1);
    }


    /* Read in the second file  */
    kname = strdup("convolvec_in2.sac");
    rsac1(kname, yarray2, &nlen2, &beg, &delta, &max, &nerr, SAC_STRING_LENGTH);

    if (nerr != 0) {
      fprintf(stderr, "Error reading in file: %s\n", kname);
      exit(-1);
    }

    /* Reverse the First Signal */
    j = 0;
    for(i = nlen1 - 1; i >= 0; i--) {
      yarray1[j] = ytmp[i];
      j++;
    }
    
    nlen = nlen1;
    if(nlen2 > nlen) {
      nlen = nlen2;
    }
    /* Allocate space for the correlation of yarray1 and yarray2 */
    max = next2((2 * nlen) - 1) * 2;
    out = (float *) malloc(sizeof(float) * max);
    if(out == NULL) {
      fprintf(stderr, "Error allocating memory for correlation\n");
      exit(-1);
    }

    /* Set up values for the cross correlation */
    nwin = 1;
    wlen = nlen;
    nfft = 0;
    
    /*     Call crscor ( Cross Correlation, no, wait, uh Convolution )
     *        - yarray1 - First  Input array to correlate
     *        - yarray2 - Second Input array to correlate
     *        - nlen    - Number of points in yarray and yarray2
     *        - nwin    - Windows to use in the correlation
     *        - wlen    - Length of the windows
     *        - type    - Type of Window (SAC_RECTANGLE)
     *        - out     - output sequence 
     *        - nfft    - Length of the output sequence
     *        - error   - Error Message
     *        - err_len - Length of Error Message (on input)
     */
    crscor(yarray1, yarray2, nlen, 
           nwin, wlen, SAC_RECTANGLE,
           out, &nfft, error, ERROR_MAX);

    /* Zero out the tmp signal */
    for(i = 0; i < MAX; i++) {
      ytmp[i] = 0.0;
    }
    
    /* Reconstruct the signal from the "cross correlation" back to front
     *  
     *  ytmp[0         : nlen1 - 2         ] <- out[nfft-nlen1+1 : nfft  - 1 ] 
     *  ytmp[nlen1 - 1 : nlen1 + nlen2  -2 ] <- out[0            : nlen2 - 1 ]
     *
     *  nfft-1 is the last point of the output sequence
     */
    for(i = 0; i <= nlen1 - 2; i++) {  
      ytmp[i] = out[nfft - nlen1 + i + 1];
    }
    for(i = 0; i <= nlen2 - 1; i++) {
      ytmp[nlen1 + i - 1] = out[i];
    }

    nfft = nlen1 + nlen2 - 1;
    xarray[0] = 0;
    leven = TRUE;
    beg = 0;
    end = beg + delta * (nfft - 1);
    newhdr();
    j = 1;
    setnhv ( "npts",   &nfft,    &nerr, SAC_STRING_LENGTH);
    setfhv ( "delta",  &delta,   &nerr, SAC_STRING_LENGTH);
    setlhv ( "leven",  &leven,   &nerr, SAC_STRING_LENGTH);
    setfhv ( "b",      &beg,     &nerr, SAC_STRING_LENGTH);
    setfhv ( "e",      &end,     &nerr, SAC_STRING_LENGTH);
    setihv ( "iftype", "itime",  &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);
    setkhv ( "kcmpnm", "Q",      &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);
    setkhv ( "kstnm",  "sta",    &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);
    setnhv ( "nwfid",  &j,       &nerr, SAC_STRING_LENGTH);
    setkhv ( "kevnm",  "FUNCGEN: TRIANGLE", &nerr, SAC_STRING_LENGTH, SAC_STRING_LENGTH);

    /*   Write out the correlation function   */
    kname = strdup("convolvec_out1.sac");    
    wsac0(kname, xarray, ytmp, &nerr, SAC_STRING_LENGTH);
    if (nerr != 0) {
      fprintf(stderr, "Error writing out file: %s\n", kname);
      exit(-1);
    }
    
    free(out);

    return 0;
}