stuplo.f 38.5 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: stuplo.f,v 1.4 2001-10-09 10:27:21 forbrig Exp $
thomas.forbriger's avatar
thomas.forbriger committed
4
c
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
c This is a simple plotting tools for seismic time series in
c SFF format
c 
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
40
c  V1.16  19/07/01   * descriptive help text for time marker
thomas.forbriger's avatar
thomas.forbriger committed
41
c  V1.17  09/10/01   * allow long postscript file names
42
43
44
45
46
47
48
49
50
c 
c======================================================================
      program stuplo
c 
c variable declaration
c --------------------
c 
c version
      character*77 version, creator
thomas.forbriger's avatar
thomas.forbriger committed
51
      parameter(version='STUPLO  V1.17  plot seismic time series (SFF format)')
52
53
54
      parameter(creator='1996 by Thomas Forbriger (IfG Stuttgart)')
c parameter definitions
      integer maxsamples, maxselect, lu, maxtraces, maxchain, maxstyle
thomas.forbriger's avatar
thomas.forbriger committed
55
      parameter(maxsamples=3 000 000)
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
      parameter(maxtraces=50)
      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
      character line*80
      integer ntrim, ntraces, i, n
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)
      character*80 filename(maxtraces), firstfree(maxtraces)
      character*10 date(maxtraces)
      character*12 time(maxtraces)
      character*5 station(maxtraces)
      character*3 channel(maxtraces)
      character*4 auxid(maxtraces)
      character*6 instype(maxtraces)
      real sectime(maxtraces), dt(maxtraces) 
      real maxval(maxtraces), average(maxtraces), minval(maxtraces)
c file reading variables
      character*80 infile
      logical moretraces
c using selections
      logical useselect
      logical selection(maxselect)
c plot caption
      character captionsel*20, captionstr*250, unitslabel*40
      character annotation*40
      integer caplen
c marker
      logical optmarker
      real xmark
c plot style
      logical optgrid, optscalex, optscaley, optinline, optabstime
      logical opttbox, ftoplab, fbotlab, optfixed, optcolor, optcenter
      logical optblack, optvertlab
      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
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
c variables related to pgplot
      character*80 device, botlabel
      real timeax(maxsamples)
c commandline
      integer maxopt, lastarg, iargc
      parameter(maxopt=23)
      character*2 optid(maxopt)
thomas.forbriger's avatar
thomas.forbriger committed
129
      character*80 optarg(maxopt)
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
      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,
     &  2h-C,2h-Y,2h-R,2h-L,2h-z,2h-m,2h-v,2h-N,2h-l,2h-h,2h-V/
      data opthasarg/.TRUE.,2*.FALSE.,3*.TRUE.,4*.FALSE.,.TRUE.,2*.FALSE.,
     &  3*.TRUE.,.FALSE.,.TRUE.,2*.FALSE.,2*.TRUE.,.FALSE./
      data optarg/3hx11,2*1h-,3hfTt,1h ,1hy,4*1h-,1h*,2*1h-,2h1.,2*2h0.,1h-,
     &  2h0.,2*1h-,8h1.,1.,1.,11h1.,1.,1.,1.,1h-/
c----------------------------------------------------------------------
c give basic information and help
      line=' '
      if (iargc().eq.1) call getarg(1, line)
c----------------------------------------------------------------------
c give help information
      if ((line(1:5).eq.'-help').or.(iargc().lt.1)) then
        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]'
        print 80,'              file [list] [nc:] ...'
        print 80,'or:    stuplo -help'
        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 *,' '
        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
173
174
175
        print *,'             ''time'' is fed directly into the PGPLOT'
        print *,'             coordinate system. It has to be given in'
        print *,'             seconds on the abscissa-scale.'
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
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
384
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
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
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
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
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
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
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
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
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
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
        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'
        print *,'             default is: fTt'
        print *,'-A comment   annotate comment to caption string'
        print *,'-i           plot caption inside box'
        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.'
        print *,' '
        print *,'styles'
        print *,'------'
        print *,' '
        print *,'-C           use different colours instead of different'
        print *,'             linestyles for marking different chains.'
        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'
        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 *,' '
        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)
      read(optarg(14), *) Yfac
      if (Yfac.lt.0.01) then
        print *,'WARNING: Y-fac is set to 1. (was too small: ',Yfac,')'
        Yfac=1.
      endif
      read(optarg(15), *) xrfac
      read(optarg(16), *) xlfac
      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)
      if (optmarker) read(optarg(18), *) xmark
      verbose=optset(19)
      optblack=optset(20)
      read(optarg(21), *) lwf_std,lwf_txt,lwf_cur
      read(optarg(22), *) chf_std,chf_xlab,chf_ylab,chf_cap
      optvertlab=optset(23)
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
        call sff_TrimLen(infile,ntrim)
        if (verbose) print 81,'  opening data file ',infile(1:ntrim)
        open(lu, file=infile, err=99, status='old')
        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,
     &                     optabstime, verbose)
              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
      ch_std=chf_std*labelheight*30.
      ch_xlab=chf_xlab*labelheight*30.
      ch_ylab=chf_ylab*yvalheight
      ch_cap=chf_cap*labelheight*30.
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)
      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
            trace=firstic(ichain)-1+ipanel
            do i=1,n
              if (captionsel(i:i).eq.'f') then
                write(line,'(a)') filename(trace)
              elseif(captionsel(i:i).eq.'t') then
                write(line,'(a)') time(trace)
                if (debug) print *,'DEBUG: time ',time(trace)
              elseif(captionsel(i:i).eq.'d') then
                write(line,'(a)') date(trace)
              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
                write(line,'(a)') station(trace)
              elseif(captionsel(i:i).eq.'i') then
                write(line,'("dt=",e12.3)') dt(trace)
              elseif(captionsel(i:i).eq.'I') then
                write(line,'(a)') instype(trace)
              elseif(captionsel(i:i).eq.'c') then
                write(line,'(a)') channel(trace)
              elseif(captionsel(i:i).eq.'a') then
                write(line,'(a)') auxid(trace)
              elseif(captionsel(i:i).eq.'A') then
                write(line,'(a)') annotation
              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
                write(line,'(a)') firstfree(trace)
              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
        do i=caplen,len(captionstr)
          captionstr(i:i)=' '
        enddo
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
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
        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 pgsvp(chf_ylab*labelheight*2.5,0.99,fracbot,fractop)
        call pgsch(ch_std)
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
            botlabel='time since midnight [sec]'
            if (opttbox) botlabel='time since midnight'
          else
            botlabel='time since first sample [sec]'
            if (opttbox) botlabel='time since first sample'
          endif
        endif
c 
c determine box style and plot box
c
        if (opttbox) then
          if (fbotlab) then
            call pgtbox('ZABCTSYHO',0.0,0,'BCTS',0.0,0)
            call pgslw(lw_txt)
            call pgsch(ch_xlab)
            call pgtbox('ZN',0.0,0,' ',0.0,0)
            call pgslw(lw_std)
            call pgsch(ch_std)
          else
            call pgtbox('ZABCTSYHO',0.0,0,'BCTS',0.0,0)
          endif
        else
          if (fbotlab) then
            call pgtbox('ABCTSYHO',0.0,0,'BCTS',0.0,0)
            call pgslw(lw_txt)
            call pgsch(ch_xlab)
            call pgtbox('N',0.0,0,' ',0.0,0)
            call pgslw(lw_std)
            call pgsch(ch_std)
          else
            call pgtbox('ABCTSYHO',0.0,0,'BCTS',0.0,0)
          endif
        endif
        call pgsch(yvalheight)
        call pgslw(lw_txt)
        call pgsch(ch_ylab)
        if (optvertlab) then
          call pgtbox(' ',0.0,0,'NV',0.0,0)
        else
          call pgtbox(' ',0.0,0,'N',0.0,0)
        endif
        call pgslw(lw_std)
        call pgsch(ch_std)
        call pgsch(labelheight*30.)
        if (optgrid) then
          call pgsls(4)
          if (opttbox) then
            call pgtbox('ZG',0.0,0,'G',0.0,0)
          else
            call pgtbox('G',0.0,0,'G',0.0,0)
          endif
          call pgsls(1)
        endif
c 
c plot a marker at specified time if desired
c 
      if (optmarker) then
        call pgsave
        call pgsch(1.0)
        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)
        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
          call pgmtxt('LV',-1.,(0.9-(chf_cap-1.)*0.025),0.0,captionstr)
        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
          if (optcolor) then
            call pgsci(istyle)
          else
            call pgsls(istyle)
          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
        call pgsls(1)
        call pgsci(1)
        call pgslw(lw_std)
c enddo ipanel
      enddo
c 
c keep alive
c
      if (optkeep) then
        cursx=0.
        cursy=0.
    3   curresult=pgcurs(cursx, cursy, curchar)  
        if ((curresult.eq.1).and.(curchar.ne.'x').and.(curchar.ne.'X'))
     &    goto 3
      endif
      call pgend
      stop 
c----------------------------------------------------------------------
c formats
   80 format(a)
   81 format(2a)
   82 format(ai4)
c stop sequences
   99 stop 'ERROR: opening data file\n'
   98 stop 'ERROR: closing data file\n'
      end
      
c======================================================================
c 
c some subroutines
c
c----------------------------------------------------------------------
c
c checkmem
c
c check whether integer and real variables use the same
c amount of bytes in memory
c
      subroutine checkmem(idata, data, maxsamples)
c declare parameters
      integer maxsamples
      real data(maxsamples)
      integer idata(maxsamples)
c declare variables
      integer i,n
c go
      n=min(maxsamples,30)
      print *,'CHECKMEM: ----- START -----'
      do i=1,n
        idata(i)=i
      enddo
      write(6,'(6(i10,1x))') (idata(i), i=1,n)
      print *,' '
      do i=1,n
        data(i)=float(idata(i))
      enddo
      write(6,'(6(f10.2,1x))') (data(i), i=1,n)
      print *,' '
      write(6,'(6(i10,1x))') (idata(i), i=1,n)
      print *,' '
      do i=1,n
        idata(i)=int(data(i))
      enddo
      write(6,'(6(i10,1x))') (idata(i), i=1,n)
      print *,'CHECKMEM: ------ END ------'
      return
      end
c 
c----------------------------------------------------------------------
c
c sffopen
c
c open sff file 
c evaluate selection
c read file header
c
      subroutine sffopen(lu, filep, useselect, selection, maxselect, debug)
c declare parameters
      integer lu, filep, maxselect
      logical useselect, selection(maxselect), debug
c declare variables
      integer ierr, i
      real sffversion
      character timestamp*13, code*20, line*80
c go
c evaluate trace selections
      call getarg(filep+1, line)
      if (line(1:2).eq.'t:') then
        filep=filep+1
        useselect=.true.
        call tf_listselect(maxselect, selection, 3, line, ierr)
        if (ierr.eq.1) then
          print *,'WARNING: selection exceeds possible range',
     &            ' from 1 to',maxselect
          print *,'         selecting only up to no.',maxselect
        elseif (ierr.eq.2) then
          print *,'WARNING: missing selection list - selecting ',
     &            'all traces'
          useselect=.false.
        elseif (ierr.ne.0) then
          print *,'WARNING: unknown error code by tf_listselect'
          print *,'         selecting all traces'
          useselect=.false.
        endif
      else
        useselect=.false.
      endif
c read file header and ignore optional blocks
      call sff_RStatus(lu,sffversion,timestamp,code,ierr)
      if (debug) print *,'DEBUG: read status'
      if (ierr.ne.0) stop 'ERROR: reading status of input file\n'
c go through header
      i=1
   10 if (code(i:i).ne.' ') then
        if (code(i:i).eq.'F') then
          call sff_SkipFree(lu, ierr)
          if (debug) print *,'DEBUG: skipped FREE block of file'
          if (ierr.ne.0) stop 'ERROR: skipping FREE block\n'
        elseif (code(i:i).eq.'S') then
          read(lu, err=99, end=98, fmt='(1x)')
          if (debug) print *,'DEBUG: skipped SOURCE line of file'
        endif
        i=i+1
        goto 10
      endif
      return
   99 stop 'ERROR: skipping SOURCE line\n'
   98 stop 'ERROR: unexpected end of file when skipping SOURCE line\n'
      end
c 
c----------------------------------------------------------------------
c 
c sffread
c
c read one sff trace an skip optional blocks
c
      subroutine 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,
     &                    optabstime, verbose)
c declare parameters
      character infile*80
      logical moretraces, debug, optabstime
      integer maxsamples, lu, trace, maxtraces
      integer nsamples(maxtraces), firstsample(maxtraces)
      character*80 filename(maxtraces), firstfree(maxtraces)
      character*10 date(maxtraces)
      character*12 time(maxtraces)
      character*5 station(maxtraces)
      character*3 channel(maxtraces)
      character*4 auxid(maxtraces)
      character*6 instype(maxtraces)
      real sectime(maxtraces), dt(maxtraces)
      real maxval(maxtraces), average(maxtraces), minval(maxtraces)
      real data(maxsamples), value, minv, maxv, timeax(maxsamples)
      integer idata(maxsamples), nsamp
      double precision avg
      logical verbose
c declare variables
      integer sample, ierr, i, ntrim
      character wid2line*132, code*20, line*80
      real ampfac, stime
c go
      call sff_TrimLen(infile,ntrim)
      filename(trace)=infile(1:ntrim)
      if (trace.eq.1) then
        firstsample(trace)=1
      else
        firstsample(trace)=firstsample(trace-1)+nsamples(trace-1)
      endif
c read trace
      moretraces=.false.
      nsamp=maxsamples-firstsample(trace)
      call sff_RData(lu, wid2line, nsamp,
     &               sectime(trace), dt(trace),
     &               idata(firstsample(trace)),
     &               ampfac, code, ierr)
      nsamples(trace)=nsamp
      if ((nsamples(trace)+firstsample(trace)-1).gt.maxsamples) 
     &  stop 'ERROR: too many samples\n'
      if (ierr.ne.0) stop 'ERROR: reading trace' 
c skip optional blocks but catch first line of FREE block 
      i=1
    1 if (code(i:i).ne.' ') then
        if (code(i:i).eq.'F') then
          read (lu, err=97, end=96, fmt='(a5)') line
          if (line.ne.'FREE ') 
     &      stop 'ERROR: not a FREE block\n'
          firstfree(trace)=' '
          read (lu, err=97, end=96, fmt='(a)') line
          if (line(1:5).ne.'FREE ') then
            call sff_TrimLen(line, ntrim)
            firstfree(trace)=line(1:ntrim)
    2       read (lu, err=97, end=96, fmt='(a)') line
            if (line(1:5).ne.'FREE ') goto 2
          endif
        elseif (code(i:i).eq.'I') then
          read(lu, err=99, end=98, fmt='(1x)')
        elseif (code(i:i).eq.'D') then
          moretraces=.true.
        endif
        i=i+1
        goto 1
      endif
c translate data from integer to real
      minv=float(idata(firstsample(trace)))*ampfac
      maxv=minv
      avg=0.d0
      if (optabstime) then
        stime=sectime(trace)
      else
        stime=0.
      endif
      if (debug) print *,'DEBUG: minv:',minv,' maxv:',maxv,' avg:',avg
      do sample=firstsample(trace),(firstsample(trace)+nsamples(trace)-1)
        value=float(idata(sample))*ampfac
For faster browsing, not all history is shown. View entire blame