stuplo.f 45.3 KB
Newer Older
1
2
c this is <stuplo.f> by Thomas Forbriger 1996
c
thomas.forbriger's avatar
thomas.forbriger committed
3
c $Id$
thomas.forbriger's avatar
thomas.forbriger committed
4
c
5
6
c Copyright 1996, 2010 by Thomas Forbriger
c
thomas.forbriger's avatar
thomas.forbriger committed
7
c This is a simple plotting tool for seismic time series in
8
c SFF format
9
10
11
12
13
14
c
c ----
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version. 
15
c 
16
17
18
19
20
21
22
23
24
25
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c GNU General Public License for more details.
c 
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c ----
c
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
c (sff is the Stuttgart File Format and is defined in sff.doc)
c 
c BE SHURE YOUR COMPUTER USES 4 BYTES OF MEMORY TO STORE
c INTEGER AND REAL VARIABLES 
c (look at the definition of idata and data in the subroutine checkmem)
c
c REVISIONS and CHANGES
c  V1.0   14/11/96   first running version
c  V1.1   15/11/96   still under major construction
c  V1.2   19/11/96   first complete release
c  V1.3   20/11/96   correct labelling of time axis
c  V1.4   26/11/96   added array traceinfile
c  V1.5   03/12/96   changed control output when reading traces
c                    and added keep option
c  V1.6   13/12/96   set correct file status on open command
c  V1.7   07/01/97   corrected timeax-calculation and timeax min/max
c  V1.8   21/01/97   correct nsamp parameter to libstuff V1.07 reading rouintes
c  V1.9   23/01/97   fixed problem with traces of zero amplitude
c  V1.10  03/02/97   * inlcuded chaining for overlay plots
c                    * included scaling factor for y- and x-size
c                    * included y-centering
c  V1.11  09/07/98   introduced time marker plotting
c  V1.12  15/07/98   * do not provide usage information in case of correct
c                      commandline parameters
c                    * introduced verbose flag
c  V1.13  10/02/99   * allow non-black colour for first set
c                    * removed old tflib calls (changed to tf_ calls)
c  V1.14  23/04/99   * linewidth and character height may be set
c                      from command line
c                    * write y-axis numeric labels perpendicular to axis
c                      if desired
c  V1.15  05/11/99   * use correct TBOX option when using full time x-labels
thomas.forbriger's avatar
thomas.forbriger committed
58
c  V1.16  19/07/01   * descriptive help text for time marker
thomas.forbriger's avatar
thomas.forbriger committed
59
c  V1.17  09/10/01   * allow long postscript file names
thomas.forbriger's avatar
thomas.forbriger committed
60
61
62
63
c  V1.18  04/04/02   new options:
c                    -T set title
c                    -X set x label
c                    -S shift time axis
thomas.forbriger's avatar
thomas.forbriger committed
64
c  V1.19             new options:
thomas.forbriger's avatar
thomas.forbriger committed
65
c                    -E force pgenv and pglab
thomas.forbriger's avatar
thomas.forbriger committed
66
67
c                    -x select time interval
c                    -y select amplitude interval
thomas.forbriger's avatar
thomas.forbriger committed
68
c  V1.19a 05/04/02   -E character scaling was incomplete
thomas.forbriger's avatar
thomas.forbriger committed
69
c  V1.19b 17/06/02   longer annotation string
thomas.forbriger's avatar
thomas.forbriger committed
70
c  V1.20  26/06/02   new option -W
71
c  V1.21  28/12/03   new skipdata subroutine
thomas.forbriger's avatar
thomas.forbriger committed
72
c  V1.21a 13/01/04   improved time scale
thomas.forbriger's avatar
thomas.forbriger committed
73
c  V1.21b 27/01/04   increased title length
thomas.forbriger's avatar
thomas.forbriger committed
74
c  V1.22  04/03/04   new option -c1
thomas.forbriger's avatar
thomas.forbriger committed
75
c  V1.23  26/09/05   new option -oS (stops line style cycling)
76
c  V1.24  29/11/05   provide new option for clean captions
thomas.forbriger's avatar
thomas.forbriger committed
77
c  V1.25  29/06/07   provide station location in caption
thomas.forbriger's avatar
thomas.forbriger committed
78
c  V1.26  02/07/07   provide option to set tick intervals
thomas.forbriger's avatar
thomas.forbriger committed
79
80
c  V1.27  18/01/08   do not use write statement for annotations
c                    just assign to character variable line
thomas.forbriger's avatar
thomas.forbriger committed
81
c  V1.28  04/12/09   use correct DIN notation for units
thomas.forbriger's avatar
thomas.forbriger committed
82
c  V1.29  13/01/11   program is prepared fro libfapidxx interface
83
84
85
86
87
88
89
90
91
c 
c======================================================================
      program stuplo
c 
c variable declaration
c --------------------
c 
c version
      character*77 version, creator
92
      parameter(version=
93
     &  'STUPLO  V1.29  plot seismic time series (SFF format)')
94
95
96
      parameter(creator='1996 by Thomas Forbriger (IfG Stuttgart)')
c parameter definitions
      integer maxsamples, maxselect, lu, maxtraces, maxchain, maxstyle
thomas.forbriger's avatar
ohj    
thomas.forbriger committed
97
      parameter(maxsamples=6 000 000)
thomas.forbriger's avatar
thomas.forbriger committed
98
      parameter(maxtraces=200)
99
100
101
102
103
104
105
106
107
108
109
110
111
      parameter(maxselect=200)
      parameter(maxchain=10)
      parameter(maxstyle=5)
      parameter(lu=10)
      real capheight, labheight, numheight
      parameter(capheight=1.5, labheight=1.5, numheight=2.)
c chaining of traces
      integer nchain, ichain, npanels, ipanel, istyle
      integer firstic(maxchain), tracesic(maxchain)
c counters
      integer trace, filep, ftrace
c variables
      logical debug, verbose, goahead
thomas.forbriger's avatar
thomas.forbriger committed
112
      character line*80
113
      integer ntrim, ntraces, i, n, ierr
114
115
116
117
118
119
120
121
122
123
c here is the one big array to contain all seismic traces
c and this big array will hold real and integer data together
c the way we do this both real and integer type variables
c must allocate 4 bytes of memory each! 
      real data(maxsamples)
      integer idata(maxsamples)
      equivalence (idata, data)
c arras to hold information on data
      integer nsamples(maxtraces), firstsample(maxtraces)
      integer traceinfile(maxtraces)
thomas.forbriger's avatar
thomas.forbriger committed
124
      character*80 filename(maxtraces), firstfree(maxtraces)
125
      character*80 informat
126
127
128
129
130
131
      character*10 date(maxtraces)
      character*12 time(maxtraces)
      character*5 station(maxtraces)
      character*3 channel(maxtraces)
      character*4 auxid(maxtraces)
      character*6 instype(maxtraces)
thomas.forbriger's avatar
thomas.forbriger committed
132
133
      character*1 loccs(maxtraces)
      real locc1(maxtraces), locc2(maxtraces)
134
135
136
      real sectime(maxtraces), dt(maxtraces) 
      real maxval(maxtraces), average(maxtraces), minval(maxtraces)
c file reading variables
thomas.forbriger's avatar
thomas.forbriger committed
137
      character*200 infile
138
139
140
141
142
143
      logical moretraces
c using selections
      logical useselect
      logical selection(maxselect)
c plot caption
      character captionsel*20, captionstr*250, unitslabel*40
thomas.forbriger's avatar
thomas.forbriger committed
144
      character annotation*100
145
146
147
148
149
150
151
      integer caplen
c marker
      logical optmarker
      real xmark
c plot style
      logical optgrid, optscalex, optscaley, optinline, optabstime
      logical opttbox, ftoplab, fbotlab, optfixed, optcolor, optcenter
thomas.forbriger's avatar
thomas.forbriger committed
152
      logical optblack, optvertlab, opttitle, optxlabel, optpgenv
thomas.forbriger's avatar
thomas.forbriger committed
153
      logical optnoxlabels, optsetxrange, optsetyrange, optwhitepaper
154
      logical optcapfromfirst, optcyclestyle, optcleancaprect
155
156
157
158
      real lwf_std,lwf_txt,lwf_cur
      real chf_std,chf_xlab,chf_ylab,chf_cap
      integer lw_std,lw_txt,lw_cur
      real ch_std,ch_xlab,ch_ylab,ch_cap
thomas.forbriger's avatar
thomas.forbriger committed
159
      character*200 partitle, parxlabel
thomas.forbriger's avatar
thomas.forbriger committed
160
      real partimeoff, parxmin, parxmax, parymin, parymax
thomas.forbriger's avatar
thomas.forbriger committed
161
162
      real majorxticks,minorxticks
      integer nminorxticks,nmajorxticks
163
164
165
166
167
168
169
170
171
172
c cursor action
      logical optkeep
      integer curresult, pgcurs
      character curchar*2
      real cursx, cursy
c plot range
      real xmaxp, xminp, ymaxp, yminp, xmaxt, xmint, ymaxt, ymint
      real nlabels, Yfac, yamp, ymid, xlfac, xrfac, xamp
      real labelheight,fractop, fracbot, vptop, vpbot, boxratio
      real yvalheight
173
174
175
176
c box coordinates
      real boxxcoor(4), boxycoor(4)
c variables to do some intermediate calculations
      real jx1,jx2,jx3,jx4,jx5
177
178
179
180
181
c variables related to pgplot
      character*80 device, botlabel
      real timeax(maxsamples)
c commandline
      integer maxopt, lastarg, iargc
182
      parameter(maxopt=37)
thomas.forbriger's avatar
thomas.forbriger committed
183
      character*3 optid(maxopt)
thomas.forbriger's avatar
thomas.forbriger committed
184
      character*200 optarg(maxopt)
185
186
187
      logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
      data optid/2h-d,2h-D,2h-g,2h-c,2h-u,2h-s,2h-i,2h-a,2h-t,2h-f,2h-A,2h-k,
thomas.forbriger's avatar
thomas.forbriger committed
188
     &  2h-C,2h-Y,2h-R,2h-L,2h-z,2h-m,2h-v,2h-N,2h-l,2h-h,2h-V,2h-X,
189
190
     &  2h-T,2h-S,2h-E,3h-nT,2h-x,2h-y,2h-W,3h-Fc,3h-oS,3h-wc,3h-xt,
     &  3h-nx,'-ty'/
191
      data opthasarg/.TRUE.,2*.FALSE.,3*.TRUE.,4*.FALSE.,.TRUE.,2*.FALSE.,
thomas.forbriger's avatar
thomas.forbriger committed
192
     &  3*.TRUE.,.FALSE.,.TRUE.,2*.FALSE.,2*.TRUE.,.FALSE.,3*.TRUE.,
193
     &  2*.FALSE.,2*.TRUE.,4*.FALSE.,3*.TRUE./
194
      data optarg/3hx11,2*1h-,3hfTt,1h ,1hy,4*1h-,1h*,2*1h-,2h1.,2*2h0.,1h-,
thomas.forbriger's avatar
thomas.forbriger committed
195
     &  2h0.,2*1h-,8h1.,1.,1.,11h1.,1.,1.,1.,1h-,2*1h-,2h0.,2*1h-,
196
     &  2*5h0.,1.,4*1h-,4h0.,0,4h0,0.,'sff'/
197
198
199
200
201
202
c----------------------------------------------------------------------
c give basic information and help
      line=' '
      if (iargc().eq.1) call getarg(1, line)
c----------------------------------------------------------------------
c give help information
203
204
205
206
      if (line(1:6).eq.'-xhelp') then
        call sff_help_details
        stop
      elseif ((line(1:5).eq.'-help').or.(iargc().lt.1)) then
207
208
209
210
211
212
213
        print 80,version
        print *,creator
        print 80,'Usage: stuplo [-d dev] [-g] [-a] [-i] [-t] [-k] [-N]'
        print 80,'              [-Y fac] [-R fac] [-L fac] [-z] [-C]'
        print 80,'              [-c options] [-s x|y|xy] [-u units]'
        print 80,'              [-f] [-A comment] [-m time] [-v] [-V]'
        print 80,'              [-l std,lab,cur] [-h std,xlab,ylab,cap]'
thomas.forbriger's avatar
thomas.forbriger committed
214
        print 80,'              [-X label] [-T title] [-S sec] [-E]'
thomas.forbriger's avatar
thomas.forbriger committed
215
        print 80,'              [-nT] [-x f,t] [-y f,t] [-W] [-Fc] [-oS]'
216
        print 80,'              [-wc] [-xt i,n] [-nx n,i] [-ty f]'
217
218
        print 80,'              file [list] [nc:] ...'
        print 80,'or:    stuplo -help'
219
        print 80,'or:    stuplo -xhelp'
220
221
222
223
224
        if (iargc().lt.1) stop 'ERROR: missing parameters\n'
        print *,' '
        print *,'STUPLO plots seismic time series that are stored'
        print *,'in SFF format (Stuttgart File Format).'
        print *,' '
225
        print *,'-ty f        select input file format (see below)'
226
227
228
229
230
231
232
233
234
235
236
237
238
239
        print *,'-d dev       This option selects the pgplot output'
        print *,'             device. Information on available ouput'
        print *,'             devices is given below.'
        print *,'-v           be verbose'
        print *,'-k           keep program and window active'
        print *,'             This is useful for opening more than one'
        print *,'             pgplot window. Stop action by pressing x'
        print *,'             or the right mouse button.'
        print *,' '
        print *,'markers'
        print *,'-------'
        print *,' '
        print *,'-g           plot grid'
        print *,'-m time      plot a marker at specified time'
thomas.forbriger's avatar
thomas.forbriger committed
240
241
242
        print *,'             ''time'' is fed directly into the PGPLOT'
        print *,'             coordinate system. It has to be given in'
        print *,'             seconds on the abscissa-scale.'
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
        print *,' '
        print *,'labels'
        print *,'------'
        print *,' '
        print *,'-t           use full time labels for x-axis (tbox style)'
        print *,'             rather than counting seconds'
        print *,'-u units     define label for y-axis'
        print *,'-V           write numeric labels along left or right'
        print *,'             viewport boundary perpendicular to axis'
        print *,'             (vertically orientated) rather than parallel'
        print *,'             to the axis.'
        print *,'-c options   define caption string'
        print *,'             the options string may consist of the'
        print *,'             following letters:'
        print *,'             f   filename'
        print *,'             t   time'
        print *,'             d   date'
        print *,'             n   nsamples'
        print *,'             i   sampling interval'
        print *,'             s   station'
        print *,'             I   instrument'
        print *,'             c   channel'
        print *,'             a   GSE2.0 auxid'
        print *,'             m   minimum and maximum value'
        print *,'             F   first line of FREE block'
        print *,'             T   number of trace'
        print *,'             A   annotation string given with -A option'
thomas.forbriger's avatar
thomas.forbriger committed
270
        print *,'             L   location of station'
271
        print *,'             default is: fTt'
272
        print *,'-wc          print caption on white box, if inside panel'
thomas.forbriger's avatar
thomas.forbriger committed
273
        print *,'-Fc          create captions only from first trace in panel'
274
275
        print *,'-A comment   annotate comment to caption string'
        print *,'-i           plot caption inside box'
thomas.forbriger's avatar
thomas.forbriger committed
276
277
278
        print *,'-T title     use this string instead of the automatically'
        print *,'             generated title'
        print *,'-X label     use this label for the time axis'
thomas.forbriger's avatar
thomas.forbriger committed
279
        print *,'-nT          use this to ommit any time scale'
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
        print *,' '
        print *,'scales and ranges'
        print *,'-----------------'
        print *,' '
        print *,'-a           use absolute time scale'
        print *,'-s x|y|xy    use same scale for all boxes'
        print *,'             x   all x-axis will have the same range'
        print *,'             y   all y-axis will have the same range'
        print *,'             xy  all x-axis and y-axis will have'
        print *,'                 the same range'
        print *,'-Y fac       The calculated maximum amplitude (rule is'
        print *,'             given by -s option) will fill the fraction'
        print *,'             fac of the plot panel (default if fac=1.).'
        print *,'-R fac       The fraction of the seismogram to cut off'
        print *,'             at the end (default fac=0.).'
        print *,'-L fac       The fraction of the seismogram to cut off'
        print *,'             at the beginning (default fac=0.).'
        print *,'-z           Set zero-line always to the center of the'
        print *,'             plot panel.'
thomas.forbriger's avatar
thomas.forbriger committed
299
        print *,'-S sec       add sec (in seconds) to the time values'
thomas.forbriger's avatar
thomas.forbriger committed
300
301
        print *,'-x f,t       set time scale from f to t (s)'
        print *,'-y f,t       set amplitude scale from f to t (s)'
thomas.forbriger's avatar
thomas.forbriger committed
302
303
304
        print *,'-xt i,n      set major tickmark interval of time axis'
        print *,'             to ''i'' and the number minor tick marks'
        print *,'             between major tick marks to ''n'' '
thomas.forbriger's avatar
thomas.forbriger committed
305
306
307
308
        print *,'-nx n,i      use ''n'' major tick marks with a minor'
        print *,'             tick marks with a minor tick mark'
        print *,'             interval of ''i'' '
        print *,'             This options overrides option -xt'
309
310
311
312
        print *,' '
        print *,'styles'
        print *,'------'
        print *,' '
thomas.forbriger's avatar
thomas.forbriger committed
313
        print *,'-W           plot black on white (also on screen)'
314
315
        print *,'-C           use different colours instead of different'
        print *,'             linestyles for marking different chains.'
thomas.forbriger's avatar
thomas.forbriger committed
316
317
        print *,'-oS          use only one style for all curves (no'
        print *,'             cycling of line style or colour).'
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
        print *,'-N           use non-black colour (or linestyle) for'
        print *,'             first dataset.'
        print *,'-l std,txt,cur'
        print *,'             The default linewidth for the following'
        print *,'             elements is scaled by the given factor:'
        print *,'             txt   text strings'
        print *,'             cur   time series curve'
        print *,'             std   all others'
        print *,'-h std,xlab,ylab,cap'
        print *,'             The default character height for the following'
        print *,'             elements is scaled by the given factor:'
        print *,'             xlab  y-axis labels'
        print *,'             ylab  x-axis labels'
        print *,'             cap   caption string'
        print *,'             std   all others'
        print *,'-f           used fixed character size'
thomas.forbriger's avatar
thomas.forbriger committed
334
335
336
337
        print *,'-E           This is a workaround for difficult cases in'
        print *,'             which the -h option does not lead to the'
        print *,'             desired result. This option forces stuplo'
        print *,'             to use a single panel and pgenv and pglab.'
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
        print *,' '
        print *,' '
        print *,'Each datafile name may be followed by a list of'
        print *,'traces. This list selects a range of traces in'
        print *,'the file which will be processed. This list may'
        print *,'contain no blank (which is the separator to the'
        print *,'next filname). The traces will always be processed'
        print *,'in the order they appear in the data file.'
        print *,' '
        print *,'Examples:'
        print *,'  t:2           will select only trace 2'
        print *,'  t:4-6,2,4     will select traces 2, 4, 5 and 6'
        print *,'  t:9,8,10,14   will select traces 8, 9, 10 and 14'
        print *,' '
        print *,'For comparison purposes it is usefull to do overlay plots.'
        print *,'This may be achieved be building severl chains of seismic'
        print *,'traces. To start a new chain in the list of files use'
        print *,'''nc:'' as a separator. Using this feature the first'
        print *,'panel will contain all first files of each chain, the'
        print *,'second one all second files, etc. Each trace will be plotted'
        print *,'with a different linestyle within each panel. Using the'
        print *,'option -C there will be a different color for each trace.'
        print *,' '
        print *,'Example for building two chains of different length with'
        print *,'the datafiles ''data1'' and ''data2'':'
        print *,' '
        print *,'    stuplo data1 t:1,2 data2 t:5 nc: data1 t:4-7'
        print *,' '
        print *,'  This is what the plotting panels will contain:'
        print *,' '
        print *,'    panel  file.trace'
        print *,'        1  data1.1 data1.4'
        print *,'        2  data1.2 data1.5'
        print *,'        3  data2.5 data1.6'
        print *,'        4  data1.7'
        print *,' '
        call pgp_showdevices
        print *,' '
        print *,'Compiled array dimensions are:'
        print *,'                 total number of samples:',maxsamples
        print *,'                 (all traces together)'
        print *,'      maximum number of traces to select:',maxselect
        print *,'        maximum number of traces to hold:',maxtraces
        print *,'       maximum number of chains to build:',maxchain
        print *,' maximum number of different plot styles:',maxstyle
        print *,' '
384
        call sff_help_formats
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
        stop
      endif
c----------------------------------------------------------------------
c set standard options
      device='x11'
      captionsel='fTt'
      unitslabel=' '
      optinline=.false.
      optscalex=.false.
      optscaley=.false.
      optabstime=.false.
      opttbox=.false.
c first read the commandline
      call tf_cmdline(1, lastarg,
     &     maxopt, optid, optarg, optset, opthasarg)
      if (optset(1)) device=optarg(1)
      debug=optset(2)
      optgrid=optset(3)
      if (optset(4)) captionsel=optarg(4)
      if (optset(5)) unitslabel=optarg(5)
      if (optset(6)) then
        if (index(optarg(6),'x').gt.0) optscalex=.true.
        if (index(optarg(6),'y').gt.0) optscaley=.true.
      endif
      optinline=optset(7)
      optabstime=optset(8)
      opttbox=optset(9)
      optfixed=optset(10)
      annotation=optarg(11)
      optkeep=optset(12)
      optcolor=optset(13)
thomas.forbriger's avatar
thomas.forbriger committed
416
      read(optarg(14), *, end=96, err=97) Yfac
417
418
419
420
      if (Yfac.lt.0.01) then
        print *,'WARNING: Y-fac is set to 1. (was too small: ',Yfac,')'
        Yfac=1.
      endif
thomas.forbriger's avatar
thomas.forbriger committed
421
422
      read(optarg(15), *, end=96, err=97) xrfac
      read(optarg(16), *, end=96, err=97) xlfac
423
424
425
426
427
428
      optcenter=optset(17)
      if ((xrfac+xlfac).ge.1.) 
     &  stop 'ERROR: check setting of -R and -L (nothing will be left)'
      if (debug) print *,'DEBUG: debug messages are switched on'
      if (debug) call checkmem(idata, data, maxsamples)
      optmarker=optset(18)
thomas.forbriger's avatar
thomas.forbriger committed
429
      if (optmarker) read(optarg(18), *, end=96, err=97) xmark
430
431
      verbose=optset(19)
      optblack=optset(20)
thomas.forbriger's avatar
thomas.forbriger committed
432
433
434
      read(optarg(21), *, end=96, err=97) lwf_std,lwf_txt,lwf_cur
      read(optarg(22), *, end=96, err=97) 
     &  chf_std,chf_xlab,chf_ylab,chf_cap
435
      optvertlab=optset(23)
thomas.forbriger's avatar
thomas.forbriger committed
436
437
438
439
440
      optxlabel=optset(24)
      parxlabel=optarg(24)
      opttitle=optset(25)
      partitle=optarg(25)
      read(optarg(26), *, end=96, err=97) partimeoff
thomas.forbriger's avatar
thomas.forbriger committed
441
442
      optpgenv=optset(27)
      optnoxlabels=optset(28)
thomas.forbriger's avatar
thomas.forbriger committed
443
444
445
446
      optsetxrange=optset(29)
      read(optarg(29), *, end=96, err=97) parxmin,parxmax
      optsetyrange=optset(30)
      read(optarg(30), *, end=96, err=97) parymin,parymax
thomas.forbriger's avatar
thomas.forbriger committed
447
      optwhitepaper=optset(31)
thomas.forbriger's avatar
thomas.forbriger committed
448
      optcapfromfirst=optset(32)
thomas.forbriger's avatar
thomas.forbriger committed
449
      optcyclestyle=(.not.optset(33))
450
      optcleancaprect=optset(34)
thomas.forbriger's avatar
thomas.forbriger committed
451
452
      read(optarg(35), *, end=96, err=97) majorxticks, nminorxticks
      read(optarg(36), *, end=96, err=97) nmajorxticks, minorxticks
453
      informat=optarg(37)
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
c      device='/krm3'
c----------------------------------------------------------------------
c initialize
      do ichain=1,maxchain
        firstic(ichain)=-1
        tracesic(ichain)=0
      enddo
c----------------------------------------------------------------------
c 
c enter main program body
c -----------------------
c
c first of all read in the time series
c
      trace=0
      nchain=1
      moretraces=.false.
      filep=lastarg+1
      goahead=.true.
      if (iargc().lt.filep) stop 'ERROR: missing data file\n'
c start file loop
    1 continue
c open new file 
        call getarg(filep, infile)
c is there a new chain to start
        if (infile(1:3).eq.'nc:') then
          nchain=nchain+1
          if (nchain.gt.maxchain) stop 'ERROR: opened too many chains'
          goto 4
        endif
484
        ntrim=index(infile,' ')-1
485
        if (verbose) print 81,'  opening data file ',infile(1:ntrim)
486
487
        call sff_select_input_format(informat, ierr)
        if (ierr.ne.0) stop 'ERROR: selecting input file format'
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
        call sffopen(lu, filep, useselect, selection, maxselect, debug)
        ftrace=0
c start trace loop
    2   continue
          ftrace=ftrace+1
          if (((.not.(useselect)).or.(selection(ftrace))).and.(goahead)) then
            if (trace.lt.maxtraces) then
              trace=trace+1
              if (firstic(nchain).lt.0) then
c start a new chain with this trace
                firstic(nchain)=trace
              endif
              tracesic(nchain)=tracesic(nchain)+1
              traceinfile(trace)=ftrace
              if (verbose) print 82,'    reading trace no.',ftrace
              call sffread(lu, infile, trace, moretraces, debug,
     &                     maxsamples, maxtraces, filename, dt, timeax,
     &                     firstsample, nsamples, firstfree, date,
     &                     time, sectime, station, channel, auxid,
     &                     instype, data, idata, maxval, average, minval,
thomas.forbriger's avatar
thomas.forbriger committed
508
509
     &                     optabstime, verbose, partimeoff,
     &                     loccs, locc1, locc2)
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
              if (debug) then
                print *,'DEBUG: dataset parameters:'
                write(6,'(" DEBUG:",6(2h >,a,1h<))')
     &            date(trace),time(trace),station(trace),
     &            channel(trace),auxid(trace),instype(trace)
                print *,'DEBUG: dt:',dt(trace),' sectime:',sectime(trace)
                print *,'DEBUG: firstfree: ',firstfree(trace)
                print *,'DEBUG: file name: ',filename(trace)
                print *,'DEBUG: first sample:',firstsample(trace),
     &                  ' last sample:',(firstsample(trace)-1+
     &                                   nsamples(trace)),
     &                  ' samples:',nsamples(trace)
                print *,'DEBUG: maxval:',maxval(trace),
     &                  ' minval:',minval(trace),
     &                  ' average:',average(trace)
              endif
            else
              print *,'WARNING: reached limit of',maxtraces,' traces'
              goahead=.false.
            endif
          else
            if (verbose) print 82,'    skipping trace no.',ftrace
            call skipdata(lu, moretraces)
          endif
          if (moretraces) goto 2
c close data file
        if (verbose) print 81,'  closing data file ',infile(1:ntrim)
        close(lu, err=98)
    4   filep=filep+1
        if (filep.le.iargc()) goto 1
      ntraces=trace
c----------------------------------------------------------------------
c maximum number of plot panels will be:
      npanels=0
      do ichain=1,nchain
        npanels=max(npanels,tracesic(ichain))
      enddo
c----------------------------------------------------------------------
c 
c all traces are read
c
c prepare data for plotting now
c
c calculate axis range
c
      xminp=timeax(firstsample(1))
      xmaxp=timeax(firstsample(1)-1+nsamples(1))
      yminp=minval(1)
      ymaxp=maxval(1)
      do trace=1,ntraces
        xminp=min(xminp,timeax(firstsample(trace)))
        xmaxp=max(xmaxp,timeax(firstsample(trace)-1+nsamples(trace)))
        yminp=min(yminp,minval(trace))
        ymaxp=max(ymaxp,maxval(trace))
      enddo
c 
c determine fractioning of surface
c 
c the way of labeling:
c
c   sx: optscalex
c   il: optinline
c   cp: caption
c   nu: numeric labels
c   la: bottom label
c
c   sx    il       cp        nu        la
c
c    x     x        -      last      last
c    -     x        -       all       all
c    x     -      all       all      last
c    -     -      all       all       all
c
c we have to check headers and scale labels
      if ((optscalex).and.(optinline)) then
c no header and only one scale label line at the bottom
        nlabels=numheight+labheight
      elseif (optinline) then
c each trace has a scale label and the last one has a label at the bottom
        nlabels=npanels*(numheight+labheight)
      elseif (optscalex) then
c each trace has a caption and numeric labels and
c the last one has a label at the bottom
        nlabels=npanels*(capheight+numheight)+labheight
      else
c each trace has caption, numeric labels and bottom label
        nlabels=npanels*(capheight+numheight+labheight)
      endif
      boxratio=40./npanels
      labelheight=1./(nlabels+npanels*boxratio)
      yvalheight=labelheight*20.
      if (optfixed) then
        labelheight=1./40.
        boxratio=((1/labelheight)-nlabels)/npanels
      endif
c no viewport up to now --> set last viewport value to top
      vpbot=0.99
c 
c----------------------------------------------------------------------
c plot style
c   linewidth
      lw_txt=lwf_txt
      lw_cur=lwf_cur
      lw_std=lwf_std
c   character height
thomas.forbriger's avatar
thomas.forbriger committed
615
616
617
618
619
620
621
622
623
624
625
      if (optpgenv) then
        ch_std=chf_std
        ch_xlab=chf_xlab
        ch_ylab=chf_ylab
        ch_cap=chf_cap
      else
        ch_std=chf_std*labelheight*30.
        ch_xlab=chf_xlab*labelheight*30.
        ch_ylab=chf_ylab*yvalheight
        ch_cap=chf_cap*labelheight*30.
      endif
626
627
628
629
630
631
632
633
634
635
636
c
c----------------------------------------------------------------------
c 
c do plot
c
      call pgp_setdevice(device,1,1)
      call pgask(.false.)
      call pgscr(3, 0.,0.,1.)
      call pgscr(4, 0.,1.,0.)
      call pgslw(lw_std)
      call pgsch(ch_std)
thomas.forbriger's avatar
thomas.forbriger committed
637
638
639
640
641
642
643
644
645
c 
c prepare background
c ------------------
c
      if (optwhitepaper) then
        call pgscr(0, 1.,1.,1.)
        call pgscr(1, 0.,0.,0.)
      endif
c
646
647
648
649
650
651
652
653
654
      do ipanel=1,npanels
c
c build caption string
c
        caplen=1
        call sff_TrimLen(captionsel,n)
        do ichain=1,nchain
c is there a trace to plot in this chain?
          if (ipanel.gt.tracesic(ichain)) goto 5
thomas.forbriger's avatar
thomas.forbriger committed
655
c or skip if only first is selcted
thomas.forbriger's avatar
thomas.forbriger committed
656
          if ((ichain.ne.1).and.(optcapfromfirst)) goto 5
657
658
659
            trace=firstic(ichain)-1+ipanel
            do i=1,n
              if (captionsel(i:i).eq.'f') then
thomas.forbriger's avatar
thomas.forbriger committed
660
                line=filename(trace)
661
              elseif(captionsel(i:i).eq.'t') then
thomas.forbriger's avatar
thomas.forbriger committed
662
                line=time(trace)
663
664
                if (debug) print *,'DEBUG: time ',time(trace)
              elseif(captionsel(i:i).eq.'d') then
thomas.forbriger's avatar
thomas.forbriger committed
665
                line=date(trace)
666
667
668
669
670
671
              elseif(captionsel(i:i).eq.'T') then
                write(line,'("(",i2,")")') traceinfile(trace)
              elseif(captionsel(i:i).eq.'n') then
                write(line,'("n=",i8)') nsamples(trace)
                if (debug) print *,'DEBUG: nsamples ',nsamples(trace)
              elseif(captionsel(i:i).eq.'s') then
thomas.forbriger's avatar
thomas.forbriger committed
672
                line=station(trace)
673
674
675
              elseif(captionsel(i:i).eq.'i') then
                write(line,'("dt=",e12.3)') dt(trace)
              elseif(captionsel(i:i).eq.'I') then
thomas.forbriger's avatar
thomas.forbriger committed
676
                line=instype(trace)
677
              elseif(captionsel(i:i).eq.'c') then
thomas.forbriger's avatar
thomas.forbriger committed
678
                line=channel(trace)
679
              elseif(captionsel(i:i).eq.'a') then
thomas.forbriger's avatar
thomas.forbriger committed
680
                line=auxid(trace)
681
              elseif(captionsel(i:i).eq.'A') then
thomas.forbriger's avatar
thomas.forbriger committed
682
                line=annotation
thomas.forbriger's avatar
thomas.forbriger committed
683
684
685
686
687
688
689
690
              elseif(captionsel(i:i).eq.'L') then
                if (loccs(trace).eq.'S') then
                  write(line,'(f8.3,"N; ",f8.3,"E")') 
     &              locc1(trace), locc2(trace)
                elseif (loccs(trace).eq.'C') then
                  write(line,'("x="f8.3,"m; y=",f8.3,"m")') 
     &              locc1(trace), locc2(trace)
                else
thomas.forbriger's avatar
thomas.forbriger committed
691
                  line='no location'
thomas.forbriger's avatar
thomas.forbriger committed
692
                endif
693
694
695
696
              elseif(captionsel(i:i).eq.'m') then
                write(line,'("min/max:",e12.3,"/",e12.3)') 
     &          minval(trace),maxval(trace)
              elseif(captionsel(i:i).eq.'F') then
thomas.forbriger's avatar
thomas.forbriger committed
697
                line=firstfree(trace)
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
              else
                print *,'WARNING: unknown caption character: ',captionsel(i:i)
                line='*'
              endif
              call sff_TrimLen(line, ntrim)
              if (debug) print *,'DEBUG: ',captionsel(i:i)
     &          ,'>',line(1:ntrim),'<', line
              if ((caplen+ntrim).lt.len(captionstr)) then
                captionstr(caplen:ntrim+caplen-1)=line(1:ntrim)
                caplen=caplen+ntrim
                if ((caplen+ntrim+3).lt.len(captionstr)) then
                  captionstr(caplen:caplen+2)='   '
                  caplen=caplen+3
                endif
              else
                print *,'WARNING: caption is truncated'
              endif
c enddo i
            enddo
    5     continue
c enddo ichain
        enddo
thomas.forbriger's avatar
thomas.forbriger committed
720
721
722
723
724
725
        if (opttitle) then
          captionstr=partitle
        else
          do i=caplen,len(captionstr)
            captionstr(i:i)=' '
          enddo
thomas.forbriger's avatar
thomas.forbriger committed
726
        endif
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
c 
c determine range of box and open plot space
c
        ichain=0
    6   continue
          ichain=ichain+1
          if (ipanel.gt.tracesic(ichain)) goto 6
        trace=firstic(ichain)-1+ipanel
        xmint=timeax(firstsample(trace))
        xmaxt=timeax(firstsample(trace)-1+nsamples(trace))
        ymint=minval(trace)
        ymaxt=maxval(trace)
        do ichain=1,nchain
c is there a trace to plot in this chain?
          if (ipanel.le.tracesic(ichain)) then
            trace=firstic(ichain)-1+ipanel
            xmint=min(xmint,timeax(firstsample(trace)))
            xmaxt=max(xmaxt,timeax(firstsample(trace)-1+nsamples(trace)))
            ymint=min(ymint,minval(trace))
            ymaxt=max(ymaxt,maxval(trace))
          endif
c enddo ichain
        enddo
c may be we have to use global scales
        if (optscalex) then
          xmint=xminp
          xmaxt=xmaxp
        endif
        if (optscaley) then
          ymint=yminp
          ymaxt=ymaxp
        endif
c modify scales by selected factors
        if (optcenter) then
          ymaxt=max(abs(ymaxt),abs(ymint))
          ymint=-max(abs(ymaxt),abs(ymint))
        endif
        yamp=(ymaxt-ymint)/Yfac
        ymid=(ymaxt+ymint)/2.
        ymint=ymid-0.5*yamp
        ymaxt=ymid+0.5*yamp
        xamp=(xmaxt-xmint)
        xmint=xmint+xamp*xlfac
        xmaxt=xmaxt-xamp*xrfac
thomas.forbriger's avatar
thomas.forbriger committed
771
772
773
774
        if (optsetxrange) then
          xmint=parxmin
          xmaxt=parxmax
        endif
775
776
777
778
779
780
781
782
783
784
785
c be aware of zero amplitudes
        if (ymaxt.eq.ymint) then
          if(ymaxt.eq.0.) then
            ymaxt=.5e-30
            ymint=-.5e-30
            print *,'NOTICE: trace ',trace,' has amplitude 0.'
          else
            ymaxt=ymint+abs(.5*ymint)
            ymint=ymint-abs(.5*ymint)
          endif
        endif
thomas.forbriger's avatar
thomas.forbriger committed
786
787
788
789
        if (optsetyrange) then
          ymint=parymin
          ymaxt=parymax
        endif
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
        if (debug) print *,'DEBUG: trace, ymaxt, ymint ',trace,ymaxt,ymint
        call pgswin(xmint, xmaxt, ymint, ymaxt)
c 
c calc values for new vp
c 
c start at botom of last box
        vptop=vpbot
c is there a caption
        fractop=vptop-(capheight*labelheight)
        if (optinline) fractop=vptop
c get bottom of box
        fracbot=fractop-boxratio*labelheight
c is there a numeric label
        vpbot=fracbot
        if (.not.((optscalex).and.(optinline))) then
          vpbot=vpbot-(numheight*labelheight)
        elseif (trace.eq.ntraces) then
          vpbot=vpbot-(numheight*labelheight)
        endif
        if (.not.(optscalex)) then
          vpbot=vpbot-(labheight*labelheight)
        elseif (trace.eq.ntraces) then
          vpbot=vpbot-(labheight*labelheight)
        endif
        call pgsch(ch_std)
thomas.forbriger's avatar
thomas.forbriger committed
815
816
817
818
819
        if (optpgenv) then
          call pgvstd
        else
          call pgsvp(chf_ylab*labelheight*2.5,0.99,fracbot,fractop)
        endif
820
821
822
823
824
825
826
827
828
829
830
831
832
c
c determine labelling style
c
        fbotlab=.true.
        ftoplab=(.not.(optinline))
        if ((optinline).and.(optscalex)) then
          fbotlab=.false.
          ftoplab=.false.
        endif
        if (trace.eq.ntraces) fbotlab=.true.
        botlabel=' '
        if ((trace.eq.ntraces).or.(.not.(optscalex))) then
          if (optabstime) then
thomas.forbriger's avatar
thomas.forbriger committed
833
            botlabel='time since midnight / sec'
834
835
            if (opttbox) botlabel='time since midnight'
          else
thomas.forbriger's avatar
thomas.forbriger committed
836
            botlabel='time since first sample / sec'
837
838
            if (opttbox) botlabel='time since first sample'
          endif
thomas.forbriger's avatar
thomas.forbriger committed
839
          if (optxlabel) botlabel=parxlabel
840
841
842
843
        endif
c 
c determine box style and plot box
c
thomas.forbriger's avatar
thomas.forbriger committed
844
845
846
847
848
849
        if (nmajorxticks.ne.0) then
          if (minorxticks.lt.((xmaxt-xmint)/1.e3))
     &      stop 'ERROR: minor tick mark interval too small!'
          nminorxticks=int((xmaxt-xmint)/(nmajorxticks*minorxticks))
          majorxticks=nminorxticks*minorxticks
        endif
850
        if (opttbox) then
thomas.forbriger's avatar
thomas.forbriger committed
851
          if (fbotlab.and.(.not.optnoxlabels)) then
thomas.forbriger's avatar
thomas.forbriger committed
852
            call pgtbox('ZABCTSHO',majorxticks,nminorxticks,'BCTS',0.0,0)
853
854
            call pgslw(lw_txt)
            call pgsch(ch_xlab)
thomas.forbriger's avatar
thomas.forbriger committed
855
            call pgtbox('ZNHO',majorxticks,nminorxticks,' ',0.0,0)
856
857
858
            call pgslw(lw_std)
            call pgsch(ch_std)
          else
thomas.forbriger's avatar
thomas.forbriger committed
859
            call pgtbox('ZABCTSHO',majorxticks,nminorxticks,'BCTS',0.0,0)
860
861
          endif
        else
thomas.forbriger's avatar
thomas.forbriger committed
862
          if (fbotlab.and.(.not.optnoxlabels)) then
thomas.forbriger's avatar
thomas.forbriger committed
863
            call pgtbox('ABCTSYHO',majorxticks,nminorxticks,'BCTS',0.0,0)
864
865
            call pgslw(lw_txt)
            call pgsch(ch_xlab)
thomas.forbriger's avatar
thomas.forbriger committed
866
            call pgtbox('NY',majorxticks,nminorxticks,' ',0.0,0)
867
868
869
            call pgslw(lw_std)
            call pgsch(ch_std)
          else
thomas.forbriger's avatar
thomas.forbriger committed
870
            call pgtbox('ABCTSYHO',majorxticks,nminorxticks,'BCTS',0.0,0)
871
872
873
874
875
876
          endif
        endif
        call pgsch(yvalheight)
        call pgslw(lw_txt)
        call pgsch(ch_ylab)
        if (optvertlab) then
thomas.forbriger's avatar
thomas.forbriger committed
877
          call pgtbox(' ',majorxticks,nminorxticks,'NV',0.0,0)
878
        else
thomas.forbriger's avatar
thomas.forbriger committed
879
          call pgtbox(' ',majorxticks,nminorxticks,'N',0.0,0)
880
881
882
883
884
885
886
        endif
        call pgslw(lw_std)
        call pgsch(ch_std)
        call pgsch(labelheight*30.)
        if (optgrid) then
          call pgsls(4)
          if (opttbox) then
thomas.forbriger's avatar
thomas.forbriger committed
887
            call pgtbox('ZG',majorxticks,nminorxticks,'G',0.0,0)
888
          else
thomas.forbriger's avatar
thomas.forbriger committed
889
            call pgtbox('G',majorxticks,nminorxticks,'G',0.0,0)
890
891
892
893
894
895
896
897
          endif
          call pgsls(1)
        endif
c 
c plot a marker at specified time if desired
c 
      if (optmarker) then
        call pgsave
thomas.forbriger's avatar
thomas.forbriger committed
898
        call pgsch(ch_std)
899
900
901
902
903
904
905
906
907
908
909
910
        call pgpt1(xmark, ymaxt, -31)
        call pgpt1(xmark, ymint, -31)
        call pgslw(3)
        call pgsls(4)
        call pgmove(xmark,ymint)
        call pgdraw(xmark,ymaxt)
        call pgunsa
      endif
c 
c plot labels
c 
        call pgslw(lw_txt)
thomas.forbriger's avatar
thomas.forbriger committed
911
912
913
914
915
916
917
918
919
        if (optpgenv) then
          call pgsch(ch_xlab)
          call pglab(botlabel, ' ', ' ')
          call pgsch(ch_ylab)
          call pglab(' ', unitslabel, ' ')
          call pgsch(ch_cap)
          if (ftoplab) then
            call pglab(' ', ' ', captionstr)
          else
920
921
922
            if (.not.optcleancaprect) then
              call pgmtxt('LV',-1.,0.9,0.0,captionstr)
            endif
thomas.forbriger's avatar
thomas.forbriger committed
923
          endif
924
        else
thomas.forbriger's avatar
thomas.forbriger committed
925
926
927
928
929
930
931
932
          call pgsch(ch_xlab)
          call pgmtxt('B',(2.9-(chf_xlab-1.)),0.5,0.5,botlabel)
          call pgsch(ch_ylab)
          call pgmtxt('L',(3.2-(chf_ylab-1.)),0.5,0.5,unitslabel)
          call pgsch(ch_cap)
          if (ftoplab) then
            call pgmtxt('T',0.7,0.5,0.5,captionstr)
          else
933
934
935
            if (.not.optcleancaprect) then
              call pgmtxt('LV',-1.,(0.9-(chf_cap-1.)*0.025),0.0,captionstr)
            endif
thomas.forbriger's avatar
thomas.forbriger committed
936
          endif
937
938
939
940
941
942
943
944
945
946
947
948
        endif
        call pgslw(lw_std)
        call pgsch(ch_std)
c 
c plot time series
c
        call pgslw(lw_cur)
        istyle=0
        if (optblack) istyle=1
        do ichain=1,nchain
          istyle=istyle+1
          if (istyle.gt.maxstyle) istyle=1
thomas.forbriger's avatar
thomas.forbriger committed
949
950
951
952
953
954
          if (optcyclestyle) then
            if (optcolor) then
              call pgsci(istyle)
            else
              call pgsls(istyle)
            endif
955
          else
thomas.forbriger's avatar
thomas.forbriger committed
956
957
958
            istyle=1
            if (optblack) istyle=2
            call pgsci(istyle)
959
960
961
962
963
964
965
966
967
968
969
          endif
c is there a trace to plot in this chain?
          if (ipanel.le.tracesic(ichain)) then
            trace=firstic(ichain)-1+ipanel
            call pgline(nsamples(trace), 
     &        timeax(firstsample(trace)),
     &        data(firstsample(trace)))
            call pgupdt
          endif
c enddo ichain 
        enddo
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
c if optcleancaprect print caption here
        if (optcleancaprect) then
c   world coordinate ranges
          jx1=xmaxt-xmint
          jx2=ymaxt-ymint
c   world coordinate offsets
          jx3=0.02*jx1
          jx4=0.02*jx2
c   first read size of box
          call pgqtxt(xmint+jx3,ymint+jx4,0.,0.,captionstr,
     &      boxxcoor, boxycoor)
          jx5=boxycoor(2)-boxycoor(1)
          call pgqtxt(xmint+jx3,ymaxt-jx4-jx5,0.,0.,captionstr,
     &      boxxcoor, boxycoor)
          call pgsave
          call pgsfs(1)
          call pgsci(0)
c   clear background of cation's box
          call pgpoly(4,boxxcoor,boxycoor)
          call pgunsa
c   write caption
          call pgptxt(xmint+jx3,ymaxt-jx4-jx5,0.,0.,captionstr)
        endif
993
994
995
996
997
998
999
1000
        call pgsls(1)
        call pgsci(1)
        call pgslw(lw_std)
c enddo ipanel
      enddo
c 
c keep alive
c
For faster browsing, not all history is shown. View entire blame