[cig-commits] r22711 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc trunk/utils/seis_process/sac2asc
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Aug 15 08:44:42 PDT 2013
Author: dkomati1
Date: 2013-08-15 08:44:42 -0700 (Thu, 15 Aug 2013)
New Revision: 22711
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac.h
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sacio.h
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac.h
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sacio.h
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/Makefile
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac2asc.c
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac2asc.c
Log:
fixed a small bug (missing subroutine argument) in sac2asc.c, added missing header files and simplified the Makefile
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/Makefile 2013-08-09 12:54:38 UTC (rev 22710)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/Makefile 2013-08-15 15:44:42 UTC (rev 22711)
@@ -1,14 +1,8 @@
#########################################################
# Makefile
-#
-# needs SAC libraries installed
########################################################
-##
-## modify to include your SAC library path
-##
-INCLUDE_SACPATH=-I/opt/seismo-util/lib/Sacio/single
-LIBRAY_SACPATH=-L/opt/seismo-util/lib
+INCLUDE_SACPATH=-I.
CFLAGS = -g -Wall $(INCLUDE_SACPATH)
@@ -16,11 +10,8 @@
OBJS=$(PROG).o
-LIBS=$(LIBRARY_SACPATH) -lSacio -lSacTools
-
-
$(PROG): $(OBJS)
- $(CC) $(CFLAGS) -o $(PROG) $(OBJS) $(LIBS)
+ $(CC) $(CFLAGS) -o $(PROG) $(OBJS)
clean:
rm -f $(PROG) $(OBJS) core *~
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac.h 2013-08-15 15:44:42 UTC (rev 22711)
@@ -0,0 +1,319 @@
+/*************************************************************************
+ Name: sac.h
+
+ Purpose: structure for header of a SAC (Seismic Analysis Code)
+ data file, and prototype for basic SAC I/O
+
+ Notes: Key to comment flags describing each field:
+ Column 1:
+ R required by SAC
+ (blank) optional
+ Column 2:
+ A = settable from a priori knowledge
+ D = available in data
+ F = available in or derivable from SEED fixed data header
+ T = available in SEED header tables
+ (blank) = not directly available from SEED data, header
+ tables, or elsewhere
+
+ Problems: none known
+
+ References: O'Neill, D. (1987). IRIS Interim Data Distribution Format
+ (SAC ASCII), Version 1.0 (12 November 1987). Incorporated
+ Research Institutions for Seismology, 1616 North Fort Myer
+ Drive, Suite 1440, Arlington, Virginia 22209. 11 pp.
+ Tull, J. (1987). SAC User's Manual, Version 10.2, October 7,
+ 1987. Lawrence Livermore National Laboratory, L-205,
+ Livermore, California 94550. ??? pp.
+
+ Language: C, hopefully ANSI standard
+
+ Author: Dennis O'Neill
+
+ Revisions: 07/15/88 Dennis O'Neill Initial preliminary release 0.9
+ 11/21/88 Dennis O'Neill Production release 1.0
+ 01/27/91 Lorraine Hwang Header number is now version 6
+ 07/06/93 Xiaoming Ding structure name sac -> sac_head
+ typedef structure to be SACHEAD
+ 12/06/96 Lupei Zhu prototype sacio functions
+**************************************************************************/
+
+
+#ifndef _sachead_h
+#define _sachead_h
+
+typedef struct sac_head
+{
+ float delta; /* RF time increment, sec */
+ float depmin; /* minimum amplitude */
+ float depmax; /* maximum amplitude */
+ float scale; /* amplitude scale factor */
+ float odelta; /* observed time inc */
+ float b; /* RD initial time - wrt nz* */
+ float e; /* RD end time */
+ float o; /* event start */
+ float a; /* 1st arrival time */
+ float internal1; /* internal use */
+ float t0; /* user-defined time pick */
+ float t1; /* user-defined time pick */
+ float t2; /* user-defined time pick */
+ float t3; /* user-defined time pick */
+ float t4; /* user-defined time pick */
+ float t5; /* user-defined time pick */
+ float t6; /* user-defined time pick */
+ float t7; /* user-defined time pick */
+ float t8; /* user-defined time pick */
+ float t9; /* user-defined time pick */
+ float f; /* event end, sec > 0 */
+ float resp0; /* instrument respnse parm*/
+ float resp1; /* instrument respnse parm*/
+ float resp2; /* instrument respnse parm*/
+ float resp3; /* instrument respnse parm*/
+ float resp4; /* instrument respnse parm*/
+ float resp5; /* instrument respnse parm*/
+ float resp6; /* instrument respnse parm*/
+ float resp7; /* instrument respnse parm*/
+ float resp8; /* instrument respnse parm*/
+ float resp9; /* instrument respnse parm*/
+ float stla; /* T station latititude */
+ float stlo; /* T station longitude */
+ float stel; /* T station elevation, m */
+ float stdp; /* T station depth, m */
+ float evla; /* event latitude */
+ float evlo; /* event longitude */
+ float evel; /* event elevation */
+ float evdp; /* event depth */
+ float unused1; /* reserved for future use*/
+ float user0; /* available to user */
+ float user1; /* available to user */
+ float user2; /* available to user */
+ float user3; /* available to user */
+ float user4; /* available to user */
+ float user5; /* available to user */
+ float user6; /* available to user */
+ float user7; /* available to user */
+ float user8; /* available to user */
+ float user9; /* available to user */
+ float dist; /* stn-event distance, km */
+ float az; /* event-stn azimuth */
+ float baz; /* stn-event azimuth */
+ float gcarc; /* stn-event dist, degrees*/
+ float internal2; /* internal use */
+ float internal3; /* internal use */
+ float depmen; /* mean value, amplitude */
+ float cmpaz; /* T component azimuth */
+ float cmpinc; /* T component inclination */
+ float unused2; /* reserved for future use*/
+ float unused3; /* reserved for future use*/
+ float unused4; /* reserved for future use*/
+ float unused5; /* reserved for future use*/
+ float unused6; /* reserved for future use*/
+ float unused7; /* reserved for future use*/
+ float unused8; /* reserved for future use*/
+ float unused9; /* reserved for future use*/
+ float unused10; /* reserved for future use*/
+ float unused11; /* reserved for future use*/
+ float unused12; /* reserved for future use*/
+ int nzyear; /* F zero time of file, yr */
+ int nzjday; /* F zero time of file, day */
+ int nzhour; /* F zero time of file, hr */
+ int nzmin; /* F zero time of file, min */
+ int nzsec; /* F zero time of file, sec */
+ int nzmsec; /* F zero time of file, msec*/
+ int internal4; /* R header version number */
+ int internal5; /* internal use */
+ int internal6; /* internal use */
+ int npts; /* RF number of samples */
+ int internal7; /* internal use */
+ int internal8; /* internal use */
+ int unused13; /* reserved for future use*/
+ int unused14; /* reserved for future use*/
+ int unused15; /* reserved for future use*/
+ int iftype; /* RA type of file */
+ int idep; /* type of amplitude */
+ int iztype; /* zero time equivalence */
+ int unused16; /* reserved for future use*/
+ int iinst; /* recording instrument */
+ int istreg; /* stn geographic region */
+ int ievreg; /* event geographic region*/
+ int ievtyp; /* event type */
+ int iqual; /* quality of data */
+ int isynth; /* synthetic data flag */
+ int unused17; /* reserved for future use*/
+ int unused18; /* reserved for future use*/
+ int unused19; /* reserved for future use*/
+ int unused20; /* reserved for future use*/
+ int unused21; /* reserved for future use*/
+ int unused22; /* reserved for future use*/
+ int unused23; /* reserved for future use*/
+ int unused24; /* reserved for future use*/
+ int unused25; /* reserved for future use*/
+ int unused26; /* reserved for future use*/
+ int leven; /* RA data-evenly-spaced flag*/
+ int lpspol; /* station polarity flag */
+ int lovrok; /* overwrite permission */
+ int lcalda; /* calc distance, azimuth */
+ int unused27; /* reserved for future use*/
+ char kstnm[8]; /* F station name */
+ char kevnm[16]; /* event name */
+ char khole[8]; /* man-made event name */
+ char ko[8]; /* event origin time id */
+ char ka[8]; /* 1st arrival time ident */
+ char kt0[8]; /* time pick 0 ident */
+ char kt1[8]; /* time pick 1 ident */
+ char kt2[8]; /* time pick 2 ident */
+ char kt3[8]; /* time pick 3 ident */
+ char kt4[8]; /* time pick 4 ident */
+ char kt5[8]; /* time pick 5 ident */
+ char kt6[8]; /* time pick 6 ident */
+ char kt7[8]; /* time pick 7 ident */
+ char kt8[8]; /* time pick 8 ident */
+ char kt9[8]; /* time pick 9 ident */
+ char kf[8]; /* end of event ident */
+ char kuser0[8]; /* available to user */
+ char kuser1[8]; /* available to user */
+ char kuser2[8]; /* available to user */
+ char kcmpnm[8]; /* F component name */
+ char knetwk[8]; /* network name */
+ char kdatrd[8]; /* date data read */
+ char kinst[8]; /* instrument name */
+} SACHEAD;
+
+
+
+/* definitions of constants for SAC enumerated data values */
+/* undocumented == couldn't find a definition for it (denio, 07/15/88) */
+#define IREAL 0 /* undocumented */
+#define ITIME 1 /* file: time series data */
+#define IRLIM 2 /* file: real&imag spectrum */
+#define IAMPH 3 /* file: ampl&phas spectrum */
+#define IXY 4 /* file: gen'l x vs y data */
+#define IUNKN 5 /* x data: unknown type */
+ /* zero time: unknown */
+ /* event type: unknown */
+#define IDISP 6 /* x data: displacement (nm) */
+#define IVEL 7 /* x data: velocity (nm/sec) */
+#define IACC 8 /* x data: accel (cm/sec/sec)*/
+#define IB 9 /* zero time: start of file */
+#define IDAY 10 /* zero time: 0000 of GMT day*/
+#define IO 11 /* zero time: event origin */
+#define IA 12 /* zero time: 1st arrival */
+#define IT0 13 /* zero time: user timepick 0*/
+#define IT1 14 /* zero time: user timepick 1*/
+#define IT2 15 /* zero time: user timepick 2*/
+#define IT3 16 /* zero time: user timepick 3*/
+#define IT4 17 /* zero time: user timepick 4*/
+#define IT5 18 /* zero time: user timepick 5*/
+#define IT6 19 /* zero time: user timepick 6*/
+#define IT7 20 /* zero time: user timepick 7*/
+#define IT8 21 /* zero time: user timepick 8*/
+#define IT9 22 /* zero time: user timepick 9*/
+#define IRADNV 23 /* undocumented */
+#define ITANNV 24 /* undocumented */
+#define IRADEV 25 /* undocumented */
+#define ITANEV 26 /* undocumented */
+#define INORTH 27 /* undocumented */
+#define IEAST 28 /* undocumented */
+#define IHORZA 29 /* undocumented */
+#define IDOWN 30 /* undocumented */
+#define IUP 31 /* undocumented */
+#define ILLLBB 32 /* undocumented */
+#define IWWSN1 33 /* undocumented */
+#define IWWSN2 34 /* undocumented */
+#define IHGLP 35 /* undocumented */
+#define ISRO 36 /* undocumented */
+#define INUCL 37 /* event type: nuclear shot */
+#define IPREN 38 /* event type: nuke pre-shot */
+#define IPOSTN 39 /* event type: nuke post-shot*/
+#define IQUAKE 40 /* event type: earthquake */
+#define IPREQ 41 /* event type: foreshock */
+#define IPOSTQ 42 /* event type: aftershock */
+#define ICHEM 43 /* event type: chemical expl */
+#define IOTHER 44 /* event type: other source */
+ /* data quality: other problm*/
+#define IGOOD 45 /* data quality: good */
+#define IGLCH 46 /* data quality: has glitches*/
+#define IDROP 47 /* data quality: has dropouts*/
+#define ILOWSN 48 /* data quality: low s/n */
+#define IRLDTA 49 /* data is real data */
+#define IVOLTS 50 /* file: velocity (volts) */
+#define INIV51 51 /* undocumented */
+#define INIV52 52 /* undocumented */
+#define INIV53 53 /* undocumented */
+#define INIV54 54 /* undocumented */
+#define INIV55 55 /* undocumented */
+#define INIV56 56 /* undocumented */
+#define INIV57 57 /* undocumented */
+#define INIV58 58 /* undocumented */
+#define INIV59 59 /* undocumented */
+#define INIV60 60 /* undocumented */
+
+/* True/false definitions */
+#ifndef TRUE
+#define FALSE 0
+#define TRUE 1
+#endif
+
+/* Format strings for writing headers for SAC ASCII files */
+#define FCS "%15.7f%15.7f%15.7f%15.7f%15.7f\n" /* for floats */
+#define ICS "%10d%10d%10d%10d%10d\n" /* for integers */
+#define CCS1 "%-8.8s%-8.8s%-8.8s\n" /* for strings */
+#define CCS2 "%-8.8s%-16.16s\n" /* for strings */
+
+
+
+/* a SAC structure containing all null values */
+static SACHEAD sac_null = {
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, 6, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }
+};
+
+/* number of bytes in header that need to be swapped on PC (int+float+long)*/
+#define HD_SIZE 440
+
+#define TMARK 10
+#define USERN 40
+
+/* prototype for SACIO functions */
+int read_sachead(const char *, SACHEAD *,int);
+float *read_sac(const char *, SACHEAD *,int);
+void sacUdata(float *, int, float *, int, int);
+int write_sac(const char *, SACHEAD, const float *,int);
+int wrtsac0(const char *, float, int, float, float, const float *,int);
+int wrtsac2(const char *, int, const float *x, const float *y,int);
+void swab4(char *, int);
+
+#endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac2asc.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac2asc.c 2013-08-09 12:54:38 UTC (rev 22710)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sac2asc.c 2013-08-15 15:44:42 UTC (rev 22711)
@@ -3,21 +3,27 @@
#include <stdio.h>
#include <sacio.h>
+#define FALSE 0
+#define TRUE 1
+
int
main(int argc, char *argv[]) {
SACHEAD sachead;
float *data;
- char* filename;
- int i;
+ char* filename;
+ int i, swap_bytes;
float time;
+/* DK DK August 2013: by default, do not swap bytes */
+ swap_bytes = FALSE;
+
if(argc < 2) {
- fprintf(stderr, "sac2asc: sacfile\n");
+ fprintf(stderr, "sac2asc: sacfile\n");
exit(-1);
}
filename = argv[1];
- if((data = read_sac(filename, &sachead) ) == 0) {
- fprintf(stderr, "Error reading sacfile: %s\n", filename);
+ if((data = read_sac(filename, &sachead, swap_bytes) ) == 0) {
+ fprintf(stderr, "Error reading sacfile: %s\n", filename);
exit(-1);
}
time = sachead.b;
@@ -29,3 +35,4 @@
return(0);
}
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sacio.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sacio.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/seis_process/sac2asc/sacio.h 2013-08-15 15:44:42 UTC (rev 22711)
@@ -0,0 +1,384 @@
+/*******************************************************************
+* sacio.c
+* SAC I/O functions:
+* read_sachead read SAC header
+* read_sac read SAC binary data
+* write_sac write SAC binary data
+* wrtsac0 write 1D array as evenly-spaced SAC binary data
+* wrtsac0_ fortran wrap for wrtsac0
+* wrtsac2 write 2 1D arrays as XY SAC data
+* wrtsac2_ fortran wrap for wrtsac2
+* swab4 reverse byte order for integer/float
+*********************************************************************/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "sac.h"
+
+/***********************************************************
+
+ read_sachead
+
+ Description: read binary SAC header from file.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD *hd SAC header to be filled
+
+ Return: 0 if success, -1 if failed
+
+ Modify history:
+ 05/29/97 Lupei Zhu Initial coding
+************************************************************/
+
+int read_sachead(const char *name,
+ SACHEAD *hd,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+
+ if ((strm = fopen(name, "rb")) == NULL) {
+ fprintf(stderr, "Unable to open %s\n",name);
+ return -1;
+ }
+
+ if (fread(hd, sizeof(SACHEAD), 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC header %s\n",name);
+ fclose(strm);
+ return -1;
+ }
+
+if(swap_bytes)
+ swab4((char *) hd, HD_SIZE);
+
+ fclose(strm);
+ return 0;
+
+}
+
+
+/***********************************************************
+
+ read_sac
+
+ Description: read binary data from file. If succeed, it will return
+ a float pointer to the read data array. The SAC header
+ is also filled. A NULL pointer is returned if failed.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD *hd SAC header to be filled
+
+ Return: float pointer to the data array, NULL if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+float* read_sac(const char *name,
+ SACHEAD *hd,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+ float *ar;
+ unsigned sz;
+
+ if ((strm = fopen(name, "rb")) == NULL) {
+ fprintf(stderr, "Unable to open %s\n",name);
+ return NULL;
+ }
+
+ if (fread(hd, sizeof(SACHEAD), 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC header %s\n",name);
+ return NULL;
+ }
+
+if(swap_bytes)
+ swab4((char *) hd, HD_SIZE);
+
+if(hd->npts < 0 || hd->npts > 1e+08 || hd->delta < 1.0e-15)
+ {
+ fprintf(stderr,"%d %e\n",hd->npts,hd->delta);
+
+ if(swap_bytes == 1)
+ swap_bytes = 0;
+ else
+ swap_bytes = 1;
+
+ swab4((char *) hd, HD_SIZE);
+ }
+
+ sz = hd->npts*sizeof(float);
+ if ((ar = (float *) malloc(sz)) == NULL) {
+ fprintf(stderr, "Error in allocating memory for reading %s\n",name);
+ return NULL;
+ }
+
+ if (fread((char *) ar, sz, 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC data %s\n",name);
+ return NULL;
+ }
+
+ fclose(strm);
+
+if(swap_bytes)
+ swab4((char *) ar, sz);
+
+ return ar;
+
+}
+
+
+
+/***********************************************************
+
+ write_sac
+
+ Description: write SAC binary data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD hd SAC header
+ const float *ar pointer to the data
+
+ Return: 0 if succeed; -1 if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+int write_sac(const char *name,
+ SACHEAD hd,
+ const float *ar,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+ unsigned sz;
+ float *data;
+ int error = 0;
+
+ sz = hd.npts*sizeof(float);
+ if (hd.iftype == IXY) sz *= 2;
+
+ if ((data = (float *) malloc(sz)) == NULL) {
+ fprintf(stderr, "Error in allocating memory for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && memcpy(data, ar, sz) == NULL) {
+ fprintf(stderr, "Error in copying data for writing %s\n",name);
+ error = 1;
+ }
+
+if(swap_bytes)
+ {
+ swab4((char *) data, sz);
+ swab4((char *) &hd, HD_SIZE);
+ }
+
+ if ( !error && (strm = fopen(name, "w")) == NULL ) {
+ fprintf(stderr,"Error in opening file for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && fwrite(&hd, sizeof(SACHEAD), 1, strm) != 1 ) {
+ fprintf(stderr,"Error in writing SAC header for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && fwrite(data, sz, 1, strm) != 1 ) {
+ fprintf(stderr,"Error in writing SAC data for writing %s\n",name);
+ error = 1;
+ }
+
+ free(data);
+ fclose(strm);
+
+ return (error==0) ? 0 : -1;
+
+}
+
+
+
+/*****************************************************
+
+ swab4
+
+ Description: reverse byte order for float/integer
+
+ Author: Lupei Zhu
+
+ Arguments: char *pt pointer to byte array
+ int n number of bytes
+
+ Return: none
+
+ Modify history:
+ 12/03/96 Lupei Zhu Initial coding
+
+************************************************************/
+
+void swab4( char *pt,
+ int n
+ )
+{
+ int i;
+ char temp;
+ for(i=0;i<n;i+=4) {
+ temp = pt[i+3];
+ pt[i+3] = pt[i];
+ pt[i] = temp;
+ temp = pt[i+2];
+ pt[i+2] = pt[i+1];
+ pt[i+1] = temp;
+ }
+}
+
+
+
+/***********************************************************
+
+ wrtsac0
+
+ Description: write 1D array into evenly spaced SAC data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ float dt sampling interval
+ int ns number of points
+ float b0 starting time
+ float dist distance range
+ const float *ar data array
+
+ Return: 0 if succeed; -1 if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swab byte-order on PC
+************************************************************/
+
+int wrtsac0(const char *name,
+ float dt,
+ int ns,
+ float b0,
+ float dist,
+ const float *ar,
+ int swap_bytes
+ )
+{
+
+ SACHEAD hd = sac_null;
+
+ hd.npts = ns;
+ hd.delta = dt;
+ hd.dist = dist;
+ hd.b = b0;
+ hd.o = 0.;
+ hd.e = b0+(hd.npts-1)*hd.delta;
+ hd.iztype = IO;
+ hd.iftype = ITIME;
+ hd.leven = TRUE;
+
+ return write_sac(name, hd, ar,swap_bytes);
+
+}
+
+
+
+/***********************************************************
+
+ wrtsac2
+
+ Description: write 2 arrays into XY SAC data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ int ns number of points
+ const float *x x data array
+ const float *y y data array
+
+ Return: 0 succeed, -1 fail
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+int wrtsac2(const char *name,
+ int n,
+ const float *x,
+ const float *y,
+ int swap_bytes
+ )
+{
+ SACHEAD hd = sac_null;
+ float *ar;
+ unsigned sz;
+ int exit_code;
+
+ hd.npts = n;
+ hd.iftype = IXY;
+ hd.leven = FALSE;
+
+ sz = n*sizeof(float);
+
+ if ( (ar = (float *) malloc(2*sz)) == NULL ) {
+ fprintf(stderr, "error in allocating memory%s\n",name);
+ return -1;
+ }
+
+ if (memcpy(ar, x, sz) == NULL) {
+ fprintf(stderr, "error in copying data %s\n",name);
+ free(ar);
+ return -1;
+ }
+ if (memcpy(ar+sz, y, sz) == NULL) {
+ fprintf(stderr, "error in copying data %s\n",name);
+ free(ar);
+ return -1;
+ }
+
+ exit_code = write_sac(name, hd, ar,swap_bytes);
+
+ free(ar);
+
+ return exit_code;
+
+}
+
+
+/* write user data t[n] into sac header starting at pos (0-9)*/
+void sacUdata(float *hd, int pos, float *t, int n, int type)
+{
+ int i;
+ hd += type;
+ for(i=pos;i<n && i<10;i++) hd[i] = t[i];
+}
+
+
+/* for fortran--write evenly-spaced data */
+void wrtsac0_(const char *name, float dt, int ns, float b0, float dist, const float *ar) {
+ wrtsac0(name, dt, ns, b0, dist, ar,0);
+}
+
+/* for fortran--write x-y data */
+void wrtsac2_(const char *name, int n, const float *x, const float *y) {
+ wrtsac2(name, n, x, y,0);
+}
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/Makefile 2013-08-09 12:54:38 UTC (rev 22710)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/Makefile 2013-08-15 15:44:42 UTC (rev 22711)
@@ -1,14 +1,8 @@
#########################################################
# Makefile
-#
-# needs SAC libraries installed
########################################################
-##
-## modify to include your SAC library path
-##
-INCLUDE_SACPATH=-I/opt/seismo-util/lib/Sacio/single
-LIBRAY_SACPATH=-L/opt/seismo-util/lib
+INCLUDE_SACPATH=-I.
CFLAGS = -g -Wall $(INCLUDE_SACPATH)
@@ -16,11 +10,8 @@
OBJS=$(PROG).o
-LIBS=$(LIBRARY_SACPATH) -lSacio -lSacTools
-
-
$(PROG): $(OBJS)
- $(CC) $(CFLAGS) -o $(PROG) $(OBJS) $(LIBS)
+ $(CC) $(CFLAGS) -o $(PROG) $(OBJS)
clean:
rm -f $(PROG) $(OBJS) core *~
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac.h 2013-08-15 15:44:42 UTC (rev 22711)
@@ -0,0 +1,319 @@
+/*************************************************************************
+ Name: sac.h
+
+ Purpose: structure for header of a SAC (Seismic Analysis Code)
+ data file, and prototype for basic SAC I/O
+
+ Notes: Key to comment flags describing each field:
+ Column 1:
+ R required by SAC
+ (blank) optional
+ Column 2:
+ A = settable from a priori knowledge
+ D = available in data
+ F = available in or derivable from SEED fixed data header
+ T = available in SEED header tables
+ (blank) = not directly available from SEED data, header
+ tables, or elsewhere
+
+ Problems: none known
+
+ References: O'Neill, D. (1987). IRIS Interim Data Distribution Format
+ (SAC ASCII), Version 1.0 (12 November 1987). Incorporated
+ Research Institutions for Seismology, 1616 North Fort Myer
+ Drive, Suite 1440, Arlington, Virginia 22209. 11 pp.
+ Tull, J. (1987). SAC User's Manual, Version 10.2, October 7,
+ 1987. Lawrence Livermore National Laboratory, L-205,
+ Livermore, California 94550. ??? pp.
+
+ Language: C, hopefully ANSI standard
+
+ Author: Dennis O'Neill
+
+ Revisions: 07/15/88 Dennis O'Neill Initial preliminary release 0.9
+ 11/21/88 Dennis O'Neill Production release 1.0
+ 01/27/91 Lorraine Hwang Header number is now version 6
+ 07/06/93 Xiaoming Ding structure name sac -> sac_head
+ typedef structure to be SACHEAD
+ 12/06/96 Lupei Zhu prototype sacio functions
+**************************************************************************/
+
+
+#ifndef _sachead_h
+#define _sachead_h
+
+typedef struct sac_head
+{
+ float delta; /* RF time increment, sec */
+ float depmin; /* minimum amplitude */
+ float depmax; /* maximum amplitude */
+ float scale; /* amplitude scale factor */
+ float odelta; /* observed time inc */
+ float b; /* RD initial time - wrt nz* */
+ float e; /* RD end time */
+ float o; /* event start */
+ float a; /* 1st arrival time */
+ float internal1; /* internal use */
+ float t0; /* user-defined time pick */
+ float t1; /* user-defined time pick */
+ float t2; /* user-defined time pick */
+ float t3; /* user-defined time pick */
+ float t4; /* user-defined time pick */
+ float t5; /* user-defined time pick */
+ float t6; /* user-defined time pick */
+ float t7; /* user-defined time pick */
+ float t8; /* user-defined time pick */
+ float t9; /* user-defined time pick */
+ float f; /* event end, sec > 0 */
+ float resp0; /* instrument respnse parm*/
+ float resp1; /* instrument respnse parm*/
+ float resp2; /* instrument respnse parm*/
+ float resp3; /* instrument respnse parm*/
+ float resp4; /* instrument respnse parm*/
+ float resp5; /* instrument respnse parm*/
+ float resp6; /* instrument respnse parm*/
+ float resp7; /* instrument respnse parm*/
+ float resp8; /* instrument respnse parm*/
+ float resp9; /* instrument respnse parm*/
+ float stla; /* T station latititude */
+ float stlo; /* T station longitude */
+ float stel; /* T station elevation, m */
+ float stdp; /* T station depth, m */
+ float evla; /* event latitude */
+ float evlo; /* event longitude */
+ float evel; /* event elevation */
+ float evdp; /* event depth */
+ float unused1; /* reserved for future use*/
+ float user0; /* available to user */
+ float user1; /* available to user */
+ float user2; /* available to user */
+ float user3; /* available to user */
+ float user4; /* available to user */
+ float user5; /* available to user */
+ float user6; /* available to user */
+ float user7; /* available to user */
+ float user8; /* available to user */
+ float user9; /* available to user */
+ float dist; /* stn-event distance, km */
+ float az; /* event-stn azimuth */
+ float baz; /* stn-event azimuth */
+ float gcarc; /* stn-event dist, degrees*/
+ float internal2; /* internal use */
+ float internal3; /* internal use */
+ float depmen; /* mean value, amplitude */
+ float cmpaz; /* T component azimuth */
+ float cmpinc; /* T component inclination */
+ float unused2; /* reserved for future use*/
+ float unused3; /* reserved for future use*/
+ float unused4; /* reserved for future use*/
+ float unused5; /* reserved for future use*/
+ float unused6; /* reserved for future use*/
+ float unused7; /* reserved for future use*/
+ float unused8; /* reserved for future use*/
+ float unused9; /* reserved for future use*/
+ float unused10; /* reserved for future use*/
+ float unused11; /* reserved for future use*/
+ float unused12; /* reserved for future use*/
+ int nzyear; /* F zero time of file, yr */
+ int nzjday; /* F zero time of file, day */
+ int nzhour; /* F zero time of file, hr */
+ int nzmin; /* F zero time of file, min */
+ int nzsec; /* F zero time of file, sec */
+ int nzmsec; /* F zero time of file, msec*/
+ int internal4; /* R header version number */
+ int internal5; /* internal use */
+ int internal6; /* internal use */
+ int npts; /* RF number of samples */
+ int internal7; /* internal use */
+ int internal8; /* internal use */
+ int unused13; /* reserved for future use*/
+ int unused14; /* reserved for future use*/
+ int unused15; /* reserved for future use*/
+ int iftype; /* RA type of file */
+ int idep; /* type of amplitude */
+ int iztype; /* zero time equivalence */
+ int unused16; /* reserved for future use*/
+ int iinst; /* recording instrument */
+ int istreg; /* stn geographic region */
+ int ievreg; /* event geographic region*/
+ int ievtyp; /* event type */
+ int iqual; /* quality of data */
+ int isynth; /* synthetic data flag */
+ int unused17; /* reserved for future use*/
+ int unused18; /* reserved for future use*/
+ int unused19; /* reserved for future use*/
+ int unused20; /* reserved for future use*/
+ int unused21; /* reserved for future use*/
+ int unused22; /* reserved for future use*/
+ int unused23; /* reserved for future use*/
+ int unused24; /* reserved for future use*/
+ int unused25; /* reserved for future use*/
+ int unused26; /* reserved for future use*/
+ int leven; /* RA data-evenly-spaced flag*/
+ int lpspol; /* station polarity flag */
+ int lovrok; /* overwrite permission */
+ int lcalda; /* calc distance, azimuth */
+ int unused27; /* reserved for future use*/
+ char kstnm[8]; /* F station name */
+ char kevnm[16]; /* event name */
+ char khole[8]; /* man-made event name */
+ char ko[8]; /* event origin time id */
+ char ka[8]; /* 1st arrival time ident */
+ char kt0[8]; /* time pick 0 ident */
+ char kt1[8]; /* time pick 1 ident */
+ char kt2[8]; /* time pick 2 ident */
+ char kt3[8]; /* time pick 3 ident */
+ char kt4[8]; /* time pick 4 ident */
+ char kt5[8]; /* time pick 5 ident */
+ char kt6[8]; /* time pick 6 ident */
+ char kt7[8]; /* time pick 7 ident */
+ char kt8[8]; /* time pick 8 ident */
+ char kt9[8]; /* time pick 9 ident */
+ char kf[8]; /* end of event ident */
+ char kuser0[8]; /* available to user */
+ char kuser1[8]; /* available to user */
+ char kuser2[8]; /* available to user */
+ char kcmpnm[8]; /* F component name */
+ char knetwk[8]; /* network name */
+ char kdatrd[8]; /* date data read */
+ char kinst[8]; /* instrument name */
+} SACHEAD;
+
+
+
+/* definitions of constants for SAC enumerated data values */
+/* undocumented == couldn't find a definition for it (denio, 07/15/88) */
+#define IREAL 0 /* undocumented */
+#define ITIME 1 /* file: time series data */
+#define IRLIM 2 /* file: real&imag spectrum */
+#define IAMPH 3 /* file: ampl&phas spectrum */
+#define IXY 4 /* file: gen'l x vs y data */
+#define IUNKN 5 /* x data: unknown type */
+ /* zero time: unknown */
+ /* event type: unknown */
+#define IDISP 6 /* x data: displacement (nm) */
+#define IVEL 7 /* x data: velocity (nm/sec) */
+#define IACC 8 /* x data: accel (cm/sec/sec)*/
+#define IB 9 /* zero time: start of file */
+#define IDAY 10 /* zero time: 0000 of GMT day*/
+#define IO 11 /* zero time: event origin */
+#define IA 12 /* zero time: 1st arrival */
+#define IT0 13 /* zero time: user timepick 0*/
+#define IT1 14 /* zero time: user timepick 1*/
+#define IT2 15 /* zero time: user timepick 2*/
+#define IT3 16 /* zero time: user timepick 3*/
+#define IT4 17 /* zero time: user timepick 4*/
+#define IT5 18 /* zero time: user timepick 5*/
+#define IT6 19 /* zero time: user timepick 6*/
+#define IT7 20 /* zero time: user timepick 7*/
+#define IT8 21 /* zero time: user timepick 8*/
+#define IT9 22 /* zero time: user timepick 9*/
+#define IRADNV 23 /* undocumented */
+#define ITANNV 24 /* undocumented */
+#define IRADEV 25 /* undocumented */
+#define ITANEV 26 /* undocumented */
+#define INORTH 27 /* undocumented */
+#define IEAST 28 /* undocumented */
+#define IHORZA 29 /* undocumented */
+#define IDOWN 30 /* undocumented */
+#define IUP 31 /* undocumented */
+#define ILLLBB 32 /* undocumented */
+#define IWWSN1 33 /* undocumented */
+#define IWWSN2 34 /* undocumented */
+#define IHGLP 35 /* undocumented */
+#define ISRO 36 /* undocumented */
+#define INUCL 37 /* event type: nuclear shot */
+#define IPREN 38 /* event type: nuke pre-shot */
+#define IPOSTN 39 /* event type: nuke post-shot*/
+#define IQUAKE 40 /* event type: earthquake */
+#define IPREQ 41 /* event type: foreshock */
+#define IPOSTQ 42 /* event type: aftershock */
+#define ICHEM 43 /* event type: chemical expl */
+#define IOTHER 44 /* event type: other source */
+ /* data quality: other problm*/
+#define IGOOD 45 /* data quality: good */
+#define IGLCH 46 /* data quality: has glitches*/
+#define IDROP 47 /* data quality: has dropouts*/
+#define ILOWSN 48 /* data quality: low s/n */
+#define IRLDTA 49 /* data is real data */
+#define IVOLTS 50 /* file: velocity (volts) */
+#define INIV51 51 /* undocumented */
+#define INIV52 52 /* undocumented */
+#define INIV53 53 /* undocumented */
+#define INIV54 54 /* undocumented */
+#define INIV55 55 /* undocumented */
+#define INIV56 56 /* undocumented */
+#define INIV57 57 /* undocumented */
+#define INIV58 58 /* undocumented */
+#define INIV59 59 /* undocumented */
+#define INIV60 60 /* undocumented */
+
+/* True/false definitions */
+#ifndef TRUE
+#define FALSE 0
+#define TRUE 1
+#endif
+
+/* Format strings for writing headers for SAC ASCII files */
+#define FCS "%15.7f%15.7f%15.7f%15.7f%15.7f\n" /* for floats */
+#define ICS "%10d%10d%10d%10d%10d\n" /* for integers */
+#define CCS1 "%-8.8s%-8.8s%-8.8s\n" /* for strings */
+#define CCS2 "%-8.8s%-16.16s\n" /* for strings */
+
+
+
+/* a SAC structure containing all null values */
+static SACHEAD sac_null = {
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345., -12345., -12345., -12345., -12345.,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, 6, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ -12345, -12345, -12345, -12345, -12345,
+ { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
+ { '-','1','2','3','4','5',' ',' ' }
+};
+
+/* number of bytes in header that need to be swapped on PC (int+float+long)*/
+#define HD_SIZE 440
+
+#define TMARK 10
+#define USERN 40
+
+/* prototype for SACIO functions */
+int read_sachead(const char *, SACHEAD *,int);
+float *read_sac(const char *, SACHEAD *,int);
+void sacUdata(float *, int, float *, int, int);
+int write_sac(const char *, SACHEAD, const float *,int);
+int wrtsac0(const char *, float, int, float, float, const float *,int);
+int wrtsac2(const char *, int, const float *x, const float *y,int);
+void swab4(char *, int);
+
+#endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac2asc.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac2asc.c 2013-08-09 12:54:38 UTC (rev 22710)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sac2asc.c 2013-08-15 15:44:42 UTC (rev 22711)
@@ -3,21 +3,27 @@
#include <stdio.h>
#include <sacio.h>
+#define FALSE 0
+#define TRUE 1
+
int
main(int argc, char *argv[]) {
SACHEAD sachead;
float *data;
- char* filename;
- int i;
+ char* filename;
+ int i, swap_bytes;
float time;
+/* DK DK August 2013: by default, do not swap bytes */
+ swap_bytes = FALSE;
+
if(argc < 2) {
- fprintf(stderr, "sac2asc: sacfile\n");
+ fprintf(stderr, "sac2asc: sacfile\n");
exit(-1);
}
filename = argv[1];
- if((data = read_sac(filename, &sachead) ) == 0) {
- fprintf(stderr, "Error reading sacfile: %s\n", filename);
+ if((data = read_sac(filename, &sachead, swap_bytes) ) == 0) {
+ fprintf(stderr, "Error reading sacfile: %s\n", filename);
exit(-1);
}
time = sachead.b;
@@ -29,3 +35,4 @@
return(0);
}
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sacio.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sacio.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/seis_process/sac2asc/sacio.h 2013-08-15 15:44:42 UTC (rev 22711)
@@ -0,0 +1,384 @@
+/*******************************************************************
+* sacio.c
+* SAC I/O functions:
+* read_sachead read SAC header
+* read_sac read SAC binary data
+* write_sac write SAC binary data
+* wrtsac0 write 1D array as evenly-spaced SAC binary data
+* wrtsac0_ fortran wrap for wrtsac0
+* wrtsac2 write 2 1D arrays as XY SAC data
+* wrtsac2_ fortran wrap for wrtsac2
+* swab4 reverse byte order for integer/float
+*********************************************************************/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "sac.h"
+
+/***********************************************************
+
+ read_sachead
+
+ Description: read binary SAC header from file.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD *hd SAC header to be filled
+
+ Return: 0 if success, -1 if failed
+
+ Modify history:
+ 05/29/97 Lupei Zhu Initial coding
+************************************************************/
+
+int read_sachead(const char *name,
+ SACHEAD *hd,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+
+ if ((strm = fopen(name, "rb")) == NULL) {
+ fprintf(stderr, "Unable to open %s\n",name);
+ return -1;
+ }
+
+ if (fread(hd, sizeof(SACHEAD), 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC header %s\n",name);
+ fclose(strm);
+ return -1;
+ }
+
+if(swap_bytes)
+ swab4((char *) hd, HD_SIZE);
+
+ fclose(strm);
+ return 0;
+
+}
+
+
+/***********************************************************
+
+ read_sac
+
+ Description: read binary data from file. If succeed, it will return
+ a float pointer to the read data array. The SAC header
+ is also filled. A NULL pointer is returned if failed.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD *hd SAC header to be filled
+
+ Return: float pointer to the data array, NULL if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+float* read_sac(const char *name,
+ SACHEAD *hd,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+ float *ar;
+ unsigned sz;
+
+ if ((strm = fopen(name, "rb")) == NULL) {
+ fprintf(stderr, "Unable to open %s\n",name);
+ return NULL;
+ }
+
+ if (fread(hd, sizeof(SACHEAD), 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC header %s\n",name);
+ return NULL;
+ }
+
+if(swap_bytes)
+ swab4((char *) hd, HD_SIZE);
+
+if(hd->npts < 0 || hd->npts > 1e+08 || hd->delta < 1.0e-15)
+ {
+ fprintf(stderr,"%d %e\n",hd->npts,hd->delta);
+
+ if(swap_bytes == 1)
+ swap_bytes = 0;
+ else
+ swap_bytes = 1;
+
+ swab4((char *) hd, HD_SIZE);
+ }
+
+ sz = hd->npts*sizeof(float);
+ if ((ar = (float *) malloc(sz)) == NULL) {
+ fprintf(stderr, "Error in allocating memory for reading %s\n",name);
+ return NULL;
+ }
+
+ if (fread((char *) ar, sz, 1, strm) != 1) {
+ fprintf(stderr, "Error in reading SAC data %s\n",name);
+ return NULL;
+ }
+
+ fclose(strm);
+
+if(swap_bytes)
+ swab4((char *) ar, sz);
+
+ return ar;
+
+}
+
+
+
+/***********************************************************
+
+ write_sac
+
+ Description: write SAC binary data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ SACHEAD hd SAC header
+ const float *ar pointer to the data
+
+ Return: 0 if succeed; -1 if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+int write_sac(const char *name,
+ SACHEAD hd,
+ const float *ar,
+ int swap_bytes
+ )
+{
+ FILE *strm;
+ unsigned sz;
+ float *data;
+ int error = 0;
+
+ sz = hd.npts*sizeof(float);
+ if (hd.iftype == IXY) sz *= 2;
+
+ if ((data = (float *) malloc(sz)) == NULL) {
+ fprintf(stderr, "Error in allocating memory for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && memcpy(data, ar, sz) == NULL) {
+ fprintf(stderr, "Error in copying data for writing %s\n",name);
+ error = 1;
+ }
+
+if(swap_bytes)
+ {
+ swab4((char *) data, sz);
+ swab4((char *) &hd, HD_SIZE);
+ }
+
+ if ( !error && (strm = fopen(name, "w")) == NULL ) {
+ fprintf(stderr,"Error in opening file for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && fwrite(&hd, sizeof(SACHEAD), 1, strm) != 1 ) {
+ fprintf(stderr,"Error in writing SAC header for writing %s\n",name);
+ error = 1;
+ }
+
+ if ( !error && fwrite(data, sz, 1, strm) != 1 ) {
+ fprintf(stderr,"Error in writing SAC data for writing %s\n",name);
+ error = 1;
+ }
+
+ free(data);
+ fclose(strm);
+
+ return (error==0) ? 0 : -1;
+
+}
+
+
+
+/*****************************************************
+
+ swab4
+
+ Description: reverse byte order for float/integer
+
+ Author: Lupei Zhu
+
+ Arguments: char *pt pointer to byte array
+ int n number of bytes
+
+ Return: none
+
+ Modify history:
+ 12/03/96 Lupei Zhu Initial coding
+
+************************************************************/
+
+void swab4( char *pt,
+ int n
+ )
+{
+ int i;
+ char temp;
+ for(i=0;i<n;i+=4) {
+ temp = pt[i+3];
+ pt[i+3] = pt[i];
+ pt[i] = temp;
+ temp = pt[i+2];
+ pt[i+2] = pt[i+1];
+ pt[i+1] = temp;
+ }
+}
+
+
+
+/***********************************************************
+
+ wrtsac0
+
+ Description: write 1D array into evenly spaced SAC data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ float dt sampling interval
+ int ns number of points
+ float b0 starting time
+ float dist distance range
+ const float *ar data array
+
+ Return: 0 if succeed; -1 if failed
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swab byte-order on PC
+************************************************************/
+
+int wrtsac0(const char *name,
+ float dt,
+ int ns,
+ float b0,
+ float dist,
+ const float *ar,
+ int swap_bytes
+ )
+{
+
+ SACHEAD hd = sac_null;
+
+ hd.npts = ns;
+ hd.delta = dt;
+ hd.dist = dist;
+ hd.b = b0;
+ hd.o = 0.;
+ hd.e = b0+(hd.npts-1)*hd.delta;
+ hd.iztype = IO;
+ hd.iftype = ITIME;
+ hd.leven = TRUE;
+
+ return write_sac(name, hd, ar,swap_bytes);
+
+}
+
+
+
+/***********************************************************
+
+ wrtsac2
+
+ Description: write 2 arrays into XY SAC data.
+
+ Author: Lupei Zhu
+
+ Arguments: const char *name file name
+ int ns number of points
+ const float *x x data array
+ const float *y y data array
+
+ Return: 0 succeed, -1 fail
+
+ Modify history:
+ 09/20/93 Lupei Zhu Initial coding
+ 12/05/96 Lupei Zhu adding error handling
+ 12/06/96 Lupei Zhu swap byte-order on PC
+************************************************************/
+
+int wrtsac2(const char *name,
+ int n,
+ const float *x,
+ const float *y,
+ int swap_bytes
+ )
+{
+ SACHEAD hd = sac_null;
+ float *ar;
+ unsigned sz;
+ int exit_code;
+
+ hd.npts = n;
+ hd.iftype = IXY;
+ hd.leven = FALSE;
+
+ sz = n*sizeof(float);
+
+ if ( (ar = (float *) malloc(2*sz)) == NULL ) {
+ fprintf(stderr, "error in allocating memory%s\n",name);
+ return -1;
+ }
+
+ if (memcpy(ar, x, sz) == NULL) {
+ fprintf(stderr, "error in copying data %s\n",name);
+ free(ar);
+ return -1;
+ }
+ if (memcpy(ar+sz, y, sz) == NULL) {
+ fprintf(stderr, "error in copying data %s\n",name);
+ free(ar);
+ return -1;
+ }
+
+ exit_code = write_sac(name, hd, ar,swap_bytes);
+
+ free(ar);
+
+ return exit_code;
+
+}
+
+
+/* write user data t[n] into sac header starting at pos (0-9)*/
+void sacUdata(float *hd, int pos, float *t, int n, int type)
+{
+ int i;
+ hd += type;
+ for(i=pos;i<n && i<10;i++) hd[i] = t[i];
+}
+
+
+/* for fortran--write evenly-spaced data */
+void wrtsac0_(const char *name, float dt, int ns, float b0, float dist, const float *ar) {
+ wrtsac0(name, dt, ns, b0, dist, ar,0);
+}
+
+/* for fortran--write x-y data */
+void wrtsac2_(const char *name, int n, const float *x, const float *y) {
+ wrtsac2(name, n, x, y,0);
+}
+
More information about the CIG-COMMITS
mailing list