output_sac.c

解吸SEED格式的源代码

			p = strtok(NULL, ","); 
  
			if (p == NULL)  
			{          
				fprintf(stderr, "Error - output_sac, unable to extract event info for sac header!\n");  
				fprintf(stderr, "Event: %s\n",  
                                        summary_file_get_event()); 
   
				break;  
			} 
        
			sach.evlo = atof(p);

			/* depth */ 
			p = strtok(NULL, ","); 
  
			if (p == NULL)  
			{          
				fprintf(stderr, "Error - output_sac, unable to extract event info for sac header!\n");  
				fprintf(stderr, "Event: %s\n",  
                                        summary_file_get_event()); 
   
				break;  
			} 
       
			sach.evdp = atof(p) * 1000;  /* convert to meters */
	
			/* convert from YY/MM/DD... to time structure */
			/* do here after all strtoks are done */
			timecvt(&ev_time, event_time);
		
		}
		else
			/* event lat/lon from hypo B071 */
			if (type71_head != NULL)
			{
				sach.evla = type71_head->latitude;
				sach.evlo = type71_head->longitude;
				sach.evdp = type71_head->depth * 1000; /* convert to meters */
				timecvt(&ev_time, type71_head->origin_time);

			}
		
		if (read_summary_flag || type71_head != NULL)
		{
			int error = 0;

			error = delaz(sach.evla, sach.evlo,
				      sach.stla, sach.stlo,
				      &sach.gcarc, &sach.az, &sach.baz);


			if (error)
			{
				fprintf(stderr, "Error - unable to compute event distance/azimuth and back azimuth for station %s; channel %s!\n", 
					current_station->station, 
					current_channel->channel);

				sach.gcarc = 0; sach.az = 0; sach.baz = 0;
			}
			else
				/* 111.1 KM per Degree - gcarc is in degrees */
				sach.dist = 111.1 * sach.gcarc;

			/* timedif gives back 10000ths, convert to seconds */
			sach.o = (float) timedif(trace_time, ev_time) / 10000;
					
		}	/* if event lat/lon/depth */

		/* mysterious LLL-set values */
		sach.internal1 = 2.0;
		sach.internal4 = 6;
		sach.internal5 = 0;
		sach.internal6 = 0;
		sach.unused27 = FALSE;

		/* not sure if these are needed, but they might be. Values might be wrong */
		sach.lpspol = FALSE;
		sach.lcalda = TRUE;

/*                 +=======================================+                 */
/*=================|  build name for and open output file  |=================*/
/*                 +=======================================+                 */

		/* 1988.023.15.34.08.2800.ANMO.SPZ is the filename for a seismogram 
		 *	from year 1988, Julian day 23, 15:34:08.2800 UT, recorded 
		 *	at station ANMO from component SPZ. 
		 */

		/* if requested from WEED, place into subdirectory */
		if (EventDirFlag)
		{
			char dirname[MAXPATHLEN];

			strcpy(dirname, get_event_dir_name());

			if (!dir_exists(dirname))
			{
				if (mkdir(dirname, 0777) == -1)
				{
					fprintf(stderr, "Unable to create event directory. Using current directory!\n");
					perror("output_sac");
					strcpy(dirname, "./");
				}
 
			}
		
			chdir(dirname);
 
		}

		sprintf (outfile_name, 
			 "%04d.%03d.%02d.%02d.%02d.%04d.%s.%s.%s.%s.%c.SAC%s",
			 data_hdr->time.year,
			 data_hdr->time.day,
			 data_hdr->time.hour,
			 data_hdr->time.minute,
			 data_hdr->time.second,
			 data_hdr->time.fracsec,
			 type10.version >= 2.3 ? data_hdr->network : "",
			 data_hdr->station,
			 data_hdr->location,
			 channel,
			 input.type,
			 (ascii) ? "_ASC" : "");

		if ((outfile = fopen (outfile_name, "w+")) == NULL)
		{
			fprintf (stderr, "\tWARNING (output_sac):  ");
			fprintf (stderr, "Output file %s is not available for writing.\n",
				 outfile_name);
			fprintf (stderr, "\tExecution continuing.\n");

			chdir(orig_dir);
		
			return;
		}

/*                 +=======================================+                 */
/*=================|        write a SAC ASCII file         |=================*/
/*                 +=======================================+                 */

		if (ascii)
		{
			/* describe the file being written */
			printf ("Writing %s.%s.%s.%s, %5d samples (ASCII),",
				type10.version >= 2.3 ? data_hdr->network : "",
				data_hdr->station,
				data_hdr->location,
				channel,
				data_hdr->nsamples);

			printf (" starting %04d,%03d %02d:%02d:%02d.%04d UT\n",
				data_hdr->time.year,
				data_hdr->time.day,
				data_hdr->time.hour,
				data_hdr->time.minute,
				data_hdr->time.second,
				data_hdr->time.fracsec);

			/* write the SAC header */
			fprintf (outfile, FCS, 
				 sach.delta, 
				 sach.depmin, 
				 sach.depmax, 
				 sach.scale, 
				 sach.odelta);
			fprintf (outfile, FCS, 
				 sach.b, 
				 sach.e, 
				 sach.o, 
				 sach.a, 
				 sach.internal1);
			fprintf (outfile, FCS, 
				 sach.t0, 
				 sach.t1, 
				 sach.t2, 
				 sach.t3, 
				 sach.t4);
			fprintf (outfile, FCS, 
				 sach.t5, 
				 sach.t6, 
				 sach.t7, 
				 sach.t8, 
				 sach.t9);
			fprintf (outfile, FCS, 
				 sach.f, 
				 sach.resp0, 
				 sach.resp1, 
				 sach.resp2, 
				 sach.resp3);
			fprintf (outfile, FCS, 
				 sach.resp4, 
				 sach.resp5, 
				 sach.resp6, 
				 sach.resp7, 
				 sach.resp8);
			fprintf (outfile, FCS, 
				 sach.resp9, 
				 sach.stla, 
				 sach.stlo, 
				 sach.stel, 
				 sach.stdp);
			fprintf (outfile, FCS, 
				 sach.evla,
				 sach.evlo,
				 sach.evel,
				 sach.evdp,
				 sach.unused1);
			fprintf (outfile, FCS, 
				 sach.user0,
				 sach.user1,
				 sach.user2,
				 sach.user3,
				 sach.user4);
			fprintf (outfile, FCS, 
				 sach.user5,
				 sach.user6,
				 sach.user7,
				 sach.user8,
				 sach.user9);
			fprintf (outfile, FCS, 
				 sach.dist,
				 sach.az,
				 sach.baz,
				 sach.gcarc,
				 sach.internal2);
			fprintf (outfile, FCS, 
				 sach.internal3,
				 sach.depmen,
				 sach.cmpaz,
				 sach.cmpinc,
				 sach.unused2);
			fprintf (outfile, FCS, 
				 sach.unused3,
				 sach.unused4,
				 sach.unused5,
				 sach.unused6,
				 sach.unused7);
			fprintf (outfile, FCS, 
				 sach.unused8,
				 sach.unused9,
				 sach.unused10,
				 sach.unused11,
				 sach.unused12);
			fprintf (outfile, ICS, 
				 sach.nzyear,
				 sach.nzjday,
				 sach.nzhour,
				 sach.nzmin,
				 sach.nzsec);
			fprintf (outfile, ICS, 
				 sach.nzmsec,
				 sach.internal4,
				 sach.internal5,
				 sach.internal6,
				 sach.npts);
			fprintf (outfile, ICS, 
				 sach.internal7,
				 sach.internal8,
				 sach.unused13,
				 sach.unused14,
				 sach.unused15);
			fprintf (outfile, ICS, 
				 sach.iftype,
				 sach.idep,
				 sach.iztype,
				 sach.unused16,
				 sach.iinst);
			fprintf (outfile, ICS, 
				 sach.istreg,
				 sach.ievreg,
				 sach.ievtyp,
				 sach.iqual,
				 sach.isynth);
			fprintf (outfile, ICS, 
				 sach.unused17,
				 sach.unused18,
				 sach.unused19,
				 sach.unused20,
				 sach.unused21);
			fprintf (outfile, ICS, 
				 sach.unused22,
				 sach.unused23,
				 sach.unused24,
				 sach.unused25,
				 sach.unused26);
			fprintf (outfile, ICS, 
				 sach.leven,
				 sach.lpspol,
				 sach.lovrok,
				 sach.lcalda,
				 sach.unused27);
			fprintf (outfile, CCS2, 
				 sach.kstnm,
				 sach.kevnm);
			fprintf (outfile, CCS1, 
				 sach.khole,
				 sach.ko,
				 sach.ka);
			fprintf (outfile, CCS1, 
				 sach.kt0,
				 sach.kt1,
				 sach.kt2);
			fprintf (outfile, CCS1, 
				 sach.kt3,
				 sach.kt4,
				 sach.kt5);
			fprintf (outfile, CCS1, 
				 sach.kt6,
				 sach.kt7,
				 sach.kt8);
			fprintf (outfile, CCS1, 
				 sach.kt9,
				 sach.kf,
				 sach.kuser0);
			fprintf (outfile, CCS1, 
				 sach.kuser1,
				 sach.kuser2,
				 sach.kcmpnm);
			fprintf (outfile, CCS1, 
				 sach.knetwk,
				 sach.kdatrd,
				 sach.kinst);

			j = k*(seis_buffer_length/data_hdr->num_mux_chan);

			/* reverse the polarity of the data if necessary */
			if (reverse_flag) {
				sac_reverse_seismic_data ( offset, j, data_hdr->nsamples, outfile_name);
			}

			/* write the SAC data */
			/* try to be somewhat efficient about the number of calls to fprint */
			for (i=0; i<data_hdr->nsamples; i+=5)
			{
				if (i+5 < data_hdr->nsamples) {
				fprintf (outfile, "%15.7g%15.7g%15.7g%15.7g%15.7g\n", 
					 seismic_data[offset+i+j],
					 seismic_data[offset+i+j+1],
					 seismic_data[offset+i+j+2],
					 seismic_data[offset+i+j+3],
					 seismic_data[offset+i+j+4]);
				}
				else {
					while (i < data_hdr->nsamples) {
						fprintf (outfile, "%15.7g", 
							 seismic_data[offset+i+j]);
						i++;

					}
					fprintf (outfile, "\n");
				}
				fflush(outfile);
			}
		}

/*                 +=======================================+                 */
/*=================|       write a SAC binary file         |=================*/
/*                 +=======================================+                 */

		if (binary)
		{
			/* describe the file being written */
			printf ("Writing %s.%s.%s.%s, %5d samples (binary),",
				type10.version >= 2.3 ? data_hdr->network : "",
				data_hdr->station,
				data_hdr->location,
				channel,
				data_hdr->nsamples);
			printf (" starting %04d,%03d %02d:%02d:%02d.%04d UT\n",
				data_hdr->time.year,
				data_hdr->time.day,
				data_hdr->time.hour,
				data_hdr->time.minute,
				data_hdr->time.second,
				data_hdr->time.fracsec);

			/* write the SAC header */
			if (fwrite (&sach, sizeof (struct sac), 1, outfile) != 1)
			{
				fprintf (stderr, "\tWARNING (output_sac):  ");
				fprintf (stderr, "failed to properly write SAC header to %s.\n",
					 outfile_name);
				fprintf (stderr, "\tExecution continuing.\n");
			}
			
			j = k*(seis_buffer_length/data_hdr->num_mux_chan);

			/* reverse the polarity of the data if necessary */
			if (reverse_flag) {
				sac_reverse_seismic_data ( offset, j, data_hdr->nsamples, outfile_name);
			}

			/* write the SAC data */
			if (fwrite (&seismic_data[offset+j], data_hdr->nsamples * sizeof (float), 
				    1, outfile) != 1)
			{
				fprintf (stderr, "\tWARNING (output_sac):  ");
				fprintf (stderr, "failed to properly write SAC data to %s.\n",
					 outfile_name);

				perror("output_sac()");

				fprintf (stderr, "\tExecution continuing.\n");
			}

		}

		fclose (outfile);

		current_channel = current_channel->next;

	}		/* the BIG k loop If multiplexed */

	chdir(orig_dir);

	return;

}


/* ------------------------------------------------------------------------- */
int sac_reverse_seismic_data (offset, mux_offset, nsamples, outfile_name) 
int offset;
int mux_offset;
int nsamples;
char *outfile_name;
{
	int i, j;
	char wrkstr[200];
	for (i=0; i<nsamples; i++)
		seismic_data[offset+i+mux_offset] = -seismic_data[offset+i+mux_offset];

	sprintf(wrkstr, "Data reversal initiated for SAC file %s", outfile_name);
	rdseed_alert_msg_out(wrkstr);
	return (nsamples);
}

/* ------------------------------------------------------------------------- */
int output_polesNzeros()



{
	void cnvt_end_time();

	struct time start;
	struct time end;

	char orig_dir[MAXPATHLEN];



	struct type52 *t_52 = current_channel;
	struct type50 *t_50 = current_station;

	current_station = type50_head;

	getcwd(orig_dir, MAXPATHLEN);