C++InterfacetoTauola
GLK.f
1 *//////////////////////////////////////////////////////////////////////////////
2 *// //
3 *// Pseudo-Class GLK //
4 *// //
5 *//////////////////////////////////////////////////////////////////////////////
6 *
7 *
8 *//////////////////////////////////////////////////////////////////////////////
9 *// ======================================================================= //
10 *// ========================== _GLK_ ==================================== //
11 *// ========== General Library of histogramming/ploting utilities ========= //
12 *// ========== It is similar but not identical to HBOOK and HPLOT ========= //
13 *// ======================================================================= //
14 *//////////////////////////////////////////////////////////////////////////////
15 *// //
16 *// Version: 1.30 //
17 *// Last correction: January 1999 //
18 *// //
19 *//////////////////////////////////////////////////////////////////////////////
20 *//
21 *// ********************************************************************
22 *// * History of the package: *
23 *// * MINI-HBOOK writen by S. Jadach, Rutherford Lab. 1976 *
24 *// * Rewritten December 1989 (S.J.) and in 1997 (S.J.) *
25 *// ********************************************************************
26 *//
27 *// Installation remarks:
28 *// (1) printing backslash character depends on F77 compilator,
29 *// user may need to modify definition of BS variable in HPLCAP
30 *//
31 *// Usage of the program:
32 *// (1) In many cases names and meanings of programs and their
33 *// parameters is the same as in original CERN libraries HBOOK
34 *// (2) Unlike to original HBOOK and HPLOT, all floating parameters
35 *// of the programs are in DOUBLE PRECISION !
36 *// (3) GLK stores histograms in DOUBLE PRECISION and always with
37 *// errors. DOUBLE PRECISION storage is essential for 10**7 events statistics!
38 *// (4) Output from GLK is a picture recorded as regular a LaTeX file
39 *// with frame and curves/histograms, it is easy to change fonts
40 *// add captions, merge plots, etc. by normal editing. Finally,
41 *// picture may be inserted in any place into LaTeX source of the
42 *// article.
43 *// (5) WARNING: two-dimensional histograms are not active!!!
44 *//
45 *//////////////////////////////////////////////////////////////////////////////
46 *// List of procedures, non-user subprograms in brackets //
47 *//////////////////////////////////////////////////////////////////////////////
48 * SUBR/FUNC 1 PAR. 2 PAR. 3 PAR. 4 PAR. 5 PAR. 6 PAR.
49 * ====================================================================
50 * (GLK_Initialize) ---- ---- ---- ---- ---- ----
51 * GLK_hi INT INT ---- ---- ---- ----
52 * GLK_hie INT INT ---- ---- ---- ----
53 * GLK_Fil1 INT DBL DBL ---- ---- ----
54 * GLK_Fil2 INT DBL DBL DBL ---- ----
55 * GLK_Book1 INT CHR*80 INT DBL DBL ----
56 * (GLK_OptOut) INT INT INT INT INT INT
57 * (L.F. GLK_Exist) INT ----- ------ ---- ---- ----
58 * GLK_Idopt INT CHR*4 ----- ---- ---- ----
59 * GLK_BookFun1 INT CHR*80 INT DBL DBL DP-FUNC
60 * GLK_Idopt INT CHR*4 ----- ---- ---- ----
61 * GLK_Book2 INT CHR*80 INT DBL DBL INT DBL DBL
62 * GLK_PrintAll --- ---- ---- ---- ---- ----
63 * GLK_SetNout INT ---- ---- ---- ---- ----
64 * GLK_Print INT ---- ---- ---- ---- ----
65 * GLK_Operat INT CHR*1 INT INT DBL DBL
66 * GLK_Hinbo1 INT CHR*8 INT DBL DBL ----
67 * GLK_Unpak INT DBL(*) CHR*(*) INT --- ----
68 * GLK_Pak INT DBL(*) ---- ---- --- ----
69 * GLK_Pake INT DBL(*) ---- ---- --- ----
70 * GLK_Range1 INT DBL DBL ---- --- ----
71 * GLK_Hinbo2 INT INT DBL DBL INT DBL DBL
72 * GLK_Ymaxim INT DBL ---- ---- --- ----
73 * GLK_Yminim INT DBL ---- ---- --- ----
74 * GLK_Reset INT CHR*(*) ---- ---- --- ----
75 * GLK_Delet INT ---- ---- ---- ---- ----
76 * (GLK_Copch) CHR*80 CHR*80 ---- ---- ---- ----
77 * (GLK_hadres) INT INT ---- ---- ---- ----
78 * GLK_Hrfile INT CHR*(*) CHR*(*) ---- ---- ----
79 * GLK_Hrout INT INT CHR*8 ---- ---- ----
80 * GLK_Hrin INT INT INT ---- ---- ----
81 * GLK_Hrend CHR*(*) ---- ---- ---- ---- ----
82 * ******************* HPLOT entries ******************
83 * GLK_PlInt INT ---- ---- ---- ---- ----
84 * GLK_PlCap INT ---- ---- ---- ---- ----
85 * GLK_PlEnd ---- ---- ---- ---- ---- ----
86 * GLK_Plot INT CHR*1 CHR*1 INT ---- ----
87 * (GLK_Plfram1) INT INT INT ---- ---- ----
88 * (GLK_SAxisX) INT DBL DBL INT DBL ----
89 * (GLK_SAxisY) INT DBL DBL INT DBL ----
90 * (GLK_PlHist) INT INT DBL DBL INT INT
91 * (GLK_PlHis2) INT INT DBL DBL INT INT
92 * (GLK_PlCirc) INT INT INT DBL DBL DBL
93 * (GLK_aprof) DBL INT DBL ---- ---- ----
94 * GLK_PlSet INT DBL ---- ---- ---- ----
95 * GLK_PlTitle INT CHR*80 ---- ---- ---- ----
96 * ******************* WMONIT entries ******************
97 * GLK_WtMon INT ???
98 * *******************************************************************
99 * END OF TABLE
100 * *******************************************************************
101 * Map of memory for single histogram
102 * ----------------------------------
103 * (1-7) Header
104 * ist +1 mark 9999999999999
105 * ist +2 mark 9d12 + id*10 + 9
106 * ist +3 iflag1 9d12 + iflag1*10 +9
107 * ist +4 iflag2 9d12 + iflag2*10 +9
108 * ist +5 scamin minimum y-scale
109 * ist +6 scamax maximum y-scale
110 * ist +7 jdlast address of the next histogram
111 * from previous history of calls (see hadres)
112 * ----------------------------------
113 * Binning size informations
114 * ----------------------------------
115 * One dimensional histogram Two dimensional histog.
116 * ------------------------- ----------------------
117 * (8-11) Binning information (8-15) Binning information
118 * ist2 = ist+7
119 * ist2 +1 NCHX ist2 +5 NCHY
120 * ist2 +2 XL ist2 +6 YL
121 * ist2 +3 XU ist2 +7 YU
122 * ist2 +4 FACTX ist2 +8 FACTY
123 *
124 * ----------------------------------
125 * All kind of sums and maxwt
126 * ----------------------------------
127 * One dimensional histogram Two dimensional histog.
128 * ------------------------- ----------------------
129 * (12-24) Under/over-flow average x (16-24)
130 * ist3 = ist+11
131 * ist3 +1 Underflow All nine combinations
132 * ist3 +2 Normal (U,N,O) x (U,N,O)
133 * ist3 +3 Overflow sum wt only (no errors)
134 * ist3 +4 U sum w**2
135 * ist3 +5 N sum w**2
136 * ist3 +6 O sum w**2
137 * ist3 +7 Sum 1
138 * ist3 +8 Sum wt*x
139 * ist3 +9 Sum wt*x*x
140 * ist3 +10 nevzer (GLK_WtMon)
141 * ist3 +11 nevove (GLK_WtMon)
142 * ist3 +12 nevacc (GLK_WtMon)
143 * ist3 +13 maxwt (GLK_WtMon)
144 * ----------------------------------
145 * Content of bins including errors
146 * ----------------------------------
147 * (25 to 24+2*nchx) (25 to 24 +nchx*nchy)
148 * sum wt and sum wt**2 sum wt only (no errors)
149 * ----------------------------------------------------------------
150 *//////////////////////////////////////////////////////////////////////////////
151 
152  SUBROUTINE glk_initialize
153 * *************************
154 * First Initialization called from may routines
155 * *************************************
156  IMPLICIT NONE
157 *----------------------------------------------------------------------
158  include 'GLK.h'
159  SAVE
160 *----------------------------------------------------------------------
161 * Note that backslash definition is varying from one
162 * instalation/compiler to another, you have to figure out by yourself
163 * how to fill backslash code into m_BS
164  CHARACTER*1 bbs1, bbs2
165  DATA bbs1 /'\\'/ ! IBM or HP with 'f77 +B '
166  DATA bbs2 /1h\ / ! HP f77 with standard options
167 *-----------------------------------------------
168  INTEGER init,i,k
169  DATA init /0/
170 *-----------------------------------------------
171  IF(init .NE. 0) RETURN
172  init=1
173 * default output unit
174  m_out=16
175  m_length=0
176 * color
177  m_keycol=0
178 * table range
179  m_keytbr=0
180  DO k=1,3
181  m_tabran(k)=1
182  ENDDO
183 *
184  DO k=1,80
185  m_color(k:k)=' '
186  ENDDO
187  m_color(1:1)='%'
188 *
189  DO i=1,m_idmax
190  DO k=1,3
191  m_index(i,k)=0
192  ENDDO
193  DO k=1,80
194  m_titlc(i)(k:k)=' '
195  ENDDO
196  ENDDO
197  DO k=1,m_lenmb
198  m_b(k)=0d0
199  ENDDO
200 *----------------------------------------------------------------------
201  m_bs = bbs1 ! IBM or HP with 'f77 +B '
202 ** m_BS = BBS2 ! HP standard options
203 *----------------------------------------------------------------------
204  END
205 
206  SUBROUTINE glk_flush
207 * ********************
208 * FLUSH memory, all histos erased!
209 * *************************************
210  IMPLICIT NONE
211  include 'GLK.h'
212  INTEGER i,k
213 *------------------------------------------------
214  CALL glk_initialize
215  m_length=0
216  DO i=1,m_idmax
217  DO k=1,3
218  m_index(i,k)=0
219  ENDDO
220  DO k=1,80
221  m_titlc(i)(k:k)=' '
222  ENDDO
223  ENDDO
224  DO k=1,m_lenmb
225  m_b(k)=0d0
226  ENDDO
227  END
228 
229  LOGICAL FUNCTION glk_exist(id)
230 * ******************************
231 * this function is true when id exists !!!!
232 * ***************************
233  IMPLICIT NONE
234  include 'GLK.h'
235  INTEGER id,lact
236 *------------------------------------------------
237  CALL glk_hadres(id,lact)
238  glk_exist = lact .NE. 0
239 *** IF(GLK_Exist) WRITE(6,*) 'GLK_Exist: does ID,lact=',id,lact
240 *** IF(.not.GLK_Exist) write(6,*) 'GLK_Exist: doesnt ID,lact=',id,lact
241  END
242 
243  DOUBLE PRECISION FUNCTION glk_hi(id,ib)
244 * **********************
245 * getting out bin content
246 * S.J. 18-Nov. 90
247 * ***********************************
248  IMPLICIT NONE
249  include 'GLK.h'
250  INTEGER id,ib
251 * locals
252  INTEGER ist,ist2,ist3,iflag2,ityphi,nch,idmem,lact
253  SAVE idmem
254  DATA idmem / -1256765/
255 *------------------------------------------------------
256  IF(id .EQ. idmem) goto 100
257  idmem=id
258 * some checks, not repeated if id the same as previously
259  CALL glk_hadres(id,lact)
260  IF(lact .EQ. 0) THEN
261  CALL glk_stop1(' GLK_hi: nonexisting histo id=',id)
262  ENDIF
263  ist = m_index(lact,2)
264  ist2 = ist+7
265  ist3 = ist+11
266 * checking if histo is of proper type
267  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
268  ityphi = mod(iflag2,10)
269  IF(ityphi .NE. 1 .AND. ityphi.NE.3) THEN
270  CALL glk_stop1(' GLK_hi: 1-dim histos only !!! id=',id)
271  ENDIF
272  100 continue
273  nch = nint(m_b(ist2+1))
274  IF(ib .EQ. 0) THEN
275 * underflow
276  glk_hi= m_b(ist3 +1)
277  ELSEIF(ib .GE. 1.and.ib .LE. nch) THEN
278 * normal bin
279  glk_hi= m_b(ist +m_buf1+ib)
280  ELSEIF(ib .EQ. nch+1) THEN
281 * overflow
282  glk_hi= m_b(ist3 +3)
283  ELSE
284 * abnormal exit
285  CALL glk_stop1(' GLK_hi: wrong binning id=',id)
286  ENDIF
287  END
288 
289  DOUBLE PRECISION FUNCTION glk_hie(id,ib)
290 * ************************
291 * getting out error of the bin
292 * s.j. 18-nov. 90
293 * ***********************************
294  IMPLICIT NONE
295  include 'GLK.h'
296 * locals
297  INTEGER ist,ist2,ist3,iflag2,ityphi,nch,lact,ib,id
298  SAVE idmem
299  INTEGER idmem
300  DATA idmem / -1256765/
301 *---------------------------------------------------------
302  IF(id .EQ. idmem) goto 100
303  idmem=id
304 * some checks, not repeated if id the same as previously
305  CALL glk_hadres(id,lact)
306  IF(lact .EQ. 0) THEN
307  CALL glk_stop1(' GLK_hie: nonexisting histo id=',id)
308  ENDIF
309  ist = m_index(lact,2)
310  ist2 = ist+7
311  ist3 = ist+11
312 * checking if histo is of proper type
313  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
314  ityphi = mod(iflag2,10)
315  IF(ityphi .NE. 1) THEN
316  CALL glk_stop1(' GLK_hie: 1-dim histos only !!! id=',id)
317  ENDIF
318  100 CONTINUE
319  nch = m_b(ist2+1)
320  IF(ib .EQ. 0) THEN
321 * underflow
322  glk_hie= dsqrt( dabs(m_b(ist3 +4)))
323  ELSEIF(ib .GE. 1.and.ib .LE. nch) THEN
324 *...normal bin, error content
325  glk_hie= dsqrt( dabs(m_b(ist+m_buf1+nch+ib)) )
326  ELSEIF(ib .EQ. nch+1) THEN
327 * overflow
328  glk_hie= dsqrt( dabs(m_b(ist3 +6)))
329  ELSE
330 * abnormal exit
331  CALL glk_stop1('+++GLK_hie: wrong binning id= ',id)
332  ENDIF
333  END
334 
335  SUBROUTINE glk_fil1(id,xx,wtx)
336 * *****************************
337 * recommended fast filling 1-dim. histogram s.j. 18 nov. 90
338 * overflow/underflow corrected by Maciek and Zbyszek
339 * ***********************************
340  IMPLICIT NONE
341  include 'GLK.h'
342  INTEGER id
343  DOUBLE PRECISION xx,wtx
344 * local
345  INTEGER ist,ist2,ist3,iflag2,ityphi,ipose1,iposx1,kposx1,kpose1,kx,nchx,lact
346  DOUBLE PRECISION x1,wt1,xl,factx,xu
347 *-------------------------------------------------------------------------
348  CALL glk_hadres(id,lact)
349 * exit for non-existig histo
350  IF(lact .EQ. 0) RETURN
351  ist = m_index(lact,2)
352  ist2 = ist+7
353  ist3 = ist+11
354 * one-dim. histo only
355  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
356  ityphi = mod(iflag2,10)
357  IF(ityphi .NE. 1) CALL glk_stop1('+++GLK_Fil1: wrong id= ',id)
358  x1= xx
359  wt1= wtx
360  m_index(lact,3)=m_index(lact,3)+1
361 * all entries
362  m_b(ist3 +7) =m_b(ist3 +7) +1
363 * for average x
364  m_b(ist3 +8) =m_b(ist3 +8) +wt1
365  m_b(ist3 +9) =m_b(ist3 +9) +wt1*x1
366 * filling coordinates
367  nchx =m_b(ist2 +1)
368  xl =m_b(ist2 +2)
369  xu =m_b(ist2 +3)
370  factx =m_b(ist2 +4)
371  IF(x1 .LT. xl) THEN
372 * underflow
373  iposx1 = ist3 +1
374  ipose1 = ist3 +4
375  kposx1 = 0
376  ELSEIF(x1 .GT. xu) THEN
377 * or overflow
378  iposx1 = ist3 +3
379  ipose1 = ist3 +6
380  kposx1 = 0
381  ELSE
382 * or any normal bin
383  iposx1 = ist3 +2
384  ipose1 = ist3 +5
385 * or given normal bin
386  kx = (x1-xl)*factx+1d0
387  kx = min( max(kx,1) ,nchx)
388  kposx1 = ist +m_buf1+kx
389  kpose1 = ist +m_buf1+nchx+kx
390  ENDIF
391  m_b(iposx1) = m_b(iposx1) +wt1
392  m_b(ipose1) = m_b(ipose1) +wt1*wt1
393  IF( kposx1 .NE. 0) m_b(kposx1) = m_b(kposx1) +wt1
394  IF( kposx1 .NE. 0) m_b(kpose1) = m_b(kpose1) +wt1*wt1
395  END !GLK_Fil1
396 
397  SUBROUTINE glk_fil1diff(id,xx,wtx,yy,wty)
398 * *****************************************
399 * Special filling routine to fill the difference f(x)-g(y)
400 * in the case when f and g are very similar x and y are close for each event.
401 * In this case coherent filling is done if x and y fall into the same bin.
402 * Note that bin width starts to be important in this method.
403 * ***********************************
404  IMPLICIT NONE
405  include 'GLK.h'
406 *
407  INTEGER id
408  DOUBLE PRECISION xx,wtx,yy,wty
409 *
410  DOUBLE PRECISION x1,x2,wt2,wt1,factx,xl,xu
411  INTEGER ist,ist2,ist3,iflag2,ityphi,kx,ke1,ie1,kx1,kx2,ke2,ix2,ie2,nchx,lact,ix1
412 *-----------------------------------------------------------------
413  CALL glk_hadres(id,lact)
414 * exit for non-existig histo
415  IF(lact .EQ. 0) RETURN
416  ist = m_index(lact,2)
417  ist2 = ist+7
418  ist3 = ist+11
419 * one-dim. histo only
420  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
421  ityphi = mod(iflag2,10)
422  IF(ityphi .NE. 1) THEN
423  CALL glk_stop1('GLK_Fil1diff: 1-dim histos only !!! id=',id)
424  ENDIF
425  x1= xx
426  x2= yy
427  wt1= wtx
428  wt2= wty
429  m_index(lact,3)=m_index(lact,3)+1
430 * all entries
431  m_b(ist3 +7) =m_b(ist3 +7) +1
432 * for average x or y not very well defined yet
433  m_b(ist3 +8) =m_b(ist3 +8) +wt1*x1 - wt2*x2
434  m_b(ist3 +9) =m_b(ist3 +9) +wt1*x1*x1 - wt2*x2*x2
435 * filling coordinates
436  nchx =m_b(ist2 +1)
437  xl =m_b(ist2 +2)
438  xu =m_b(ist2 +3)
439  factx =m_b(ist2 +4)
440 * first variable
441  IF(x1 .LT. xl) THEN ! underflow
442  ix1 = ist3 +1
443  ie1 = ist3 +4
444  kx1 = 0
445  ELSEIF(x1 .GT. xu) THEN ! or overflow
446  ix1 = ist3 +3
447  ie1 = ist3 +6
448  kx1 = 0
449  ELSE ! normal bin
450  ix1 = ist3 +2
451  ie1 = ist3 +5
452  kx = (x1-xl)*factx+1d0
453  kx = min( max(kx,1) ,nchx)
454  kx1 = ist +m_buf1+kx
455  ke1 = ist +m_buf1+nchx+kx
456  ENDIF
457 * second variable
458  IF(x2 .LT. xl) THEN ! underflow
459  ix2 = ist3 +1
460  ie2 = ist3 +4
461  kx2 = 0
462  ELSEIF(x2 .GT. xu) THEN ! or overflow
463  ix2 = ist3 +3
464  ie2 = ist3 +6
465  kx2 = 0
466  ELSE ! normal bin
467  ix2 = ist3 +2
468  ie2 = ist3 +5
469  kx = (x2-xl)*factx+1d0
470  kx = min( max(kx,1) ,nchx)
471  kx2 = ist +m_buf1+kx
472  ke2 = ist +m_buf1+nchx+kx
473  ENDIF
474 * coherent filling
475  IF( ix1 .EQ. ix2 ) THEN
476  m_b(ix1) = m_b(ix1) +wt1-wt2
477  m_b(ie1) = m_b(ie1) +(wt1-wt2)**2
478  ELSE
479  m_b(ix1) = m_b(ix1) +wt1
480  m_b(ie1) = m_b(ie1) +wt1*wt1
481  m_b(ix2) = m_b(ix2) -wt2
482  m_b(ie2) = m_b(ie2) +wt2*wt2
483  ENDIF
484  IF( kx1 .EQ. kx2 ) THEN
485  IF( kx1 .NE. 0) THEN
486  m_b(kx1) = m_b(kx1) +wt1-wt2
487  m_b(ke1) = m_b(ke1) +(wt1-wt2)**2
488  ENDIF
489  ELSE
490  IF( kx1 .NE. 0) THEN
491  m_b(kx1) = m_b(kx1) +wt1
492  m_b(ke1) = m_b(ke1) +wt1*wt1
493  ENDIF
494  IF( kx2 .NE. 0) THEN
495  m_b(kx2) = m_b(kx2) -wt2
496  m_b(ke2) = m_b(ke2) +wt2*wt2
497  ENDIF
498  ENDIF
499  END !GLK_Fil1diff
500 
501  SUBROUTINE glk_fil2(id,x,y,wtw)
502 * ****************************
503 * this routine not finished, 1-dim only!
504 * ***********************************
505  IMPLICIT NONE
506  include 'GLK.h'
507  INTEGER id
508  DOUBLE PRECISION x,y,wtw
509 * local
510  INTEGER ist,iflag2,ityphi,ist2,ist3,nchx,nchy,ly,ky,k2,kx,lact,lx,k,l
511  DOUBLE PRECISION xx,yy,wt,factx,xl,yl,facty
512 *-------------------------------------------------------
513  CALL glk_hadres(id,lact)
514  IF(lact .EQ. 0) RETURN
515  ist = m_index(lact,2)
516 * one-dim. histo
517  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
518  ityphi = mod(iflag2,10)
519  IF(ityphi .NE. 2) THEN
520  CALL glk_stop1('GLK_Fil2: 2-dim histos only !!! id=',id)
521  ENDIF
522 *...two-dim. scattergram, no errors!
523  ist2 = ist+7
524  ist3 = ist+15
525  xx= x
526  yy= y
527  wt= wtw
528  m_index(lact,3)=m_index(lact,3)+1
529 * x-axis
530  nchx =m_b(ist2 +1)
531  xl =m_b(ist2 +2)
532  factx =m_b(ist2 +4)
533  kx=(xx-xl)*factx+1d0
534  lx=2
535  IF(kx .LT. 1) lx=1
536  IF(kx .GT. nchx) lx=3
537  l = ist+34 +lx
538  m_b(l) = m_b(l) +wt
539  k = ist+m_buf2 +kx
540  IF(lx .EQ. 2) m_b(k) =m_b(k) +wt
541  k2 = ist+m_buf2 +nchx+kx
542  IF(lx .EQ. 2) m_b(k2) =m_b(k2) +wt**2
543 * y-axix
544  nchy =m_b(ist2 +5)
545  yl =m_b(ist2 +6)
546  facty =m_b(ist2 +8)
547  ky=(yy-yl)*facty+1d0
548  ly=2
549  IF(ky .LT. 1) ly=1
550  IF(ky .GT. nchy) ly=3
551 * under/over-flow
552  l = ist3 +lx +3*(ly-1)
553  m_b(l) =m_b(l)+wt
554 * regular bin
555  k = ist+m_buf2 +kx +nchx*(ky-1)
556  IF(lx .EQ. 2.and.ly .EQ. 2) m_b(k)=m_b(k)+wt
557  END
558 
559  SUBROUTINE glk_book1(id,title,nnchx,xxl,xxu)
560 * ********************************************
561  IMPLICIT NONE
562  include 'GLK.h'
563  INTEGER id
564  DOUBLE PRECISION xxl,xxu
565  CHARACTER*80 title
566 * locals
567  DOUBLE PRECISION xl,xu,ddx
568  INTEGER ist,nchx,ioplog,iopsla,ioperb,iflag2,ityphi,iflag1
569  INTEGER ist3,ist2,lengt2,lact,nnchx,iopsc2,iopsc1,j
570  LOGICAL glk_exist
571 *-------------------------------------------------
572  CALL glk_initialize
573  IF(glk_exist(id)) goto 900
574  ist=m_length
575  CALL glk_hadres(0,lact)
576 * Check if there is a free entry in the m_index
577  IF(lact .EQ. 0)
578  $ CALL glk_stop1('GLK_Book1: to many histos !!!!!, id= ',id)
579  m_index(lact,1)=id
580  m_index(lact,2)=m_length
581  m_index(lact,3)=0
582 * -------
583  CALL glk_copch(title,m_titlc(lact))
584  nchx =nnchx
585  IF(nchx .GT. m_maxnb)
586  $ CALL glk_stop1(' GLK_Book1: To many bins requested,id= ',id)
587  xl =xxl
588  xu =xxu
589 * ---------- title and bin content ----------
590  lengt2 = m_length +2*nchx +m_buf1+1
591  IF(lengt2 .GE. m_lenmb)
592  $ CALL glk_stop1('GLK_Book1:too litle storage, m_LenmB= ',m_lenmb)
593 *
594  DO j=m_length+1,lengt2+1
595  m_b(j) = 0d0
596  ENDDO
597  m_length=lengt2
598 *... default flags
599  ioplog = 1
600  iopsla = 1
601  ioperb = 1
602  iopsc1 = 1
603  iopsc2 = 1
604  iflag1 =
605  $ ioplog+10*iopsla+100*ioperb+1000*iopsc1+10000*iopsc2
606  ityphi = 1
607  iflag2 = ityphi
608 * examples of decoding flags
609 * id = nint(m_b(ist+2)-9d0-9d12)/10
610 * iflag1 = nint(m_b(ist+3)-9d0-9d12)/10
611 * ioplog = mod(iflag1,10)
612 * iopsla = mod(iflag1,100)/10
613 * ioperb = mod(iflag1,1000)/100
614 * iopsc1 = mod(iflag1,10000)/1000
615 * iopsc2 = mod(iflag1,100000)/10000
616 * iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
617 * ityphi = mod(iflag2,10)
618 *--------- buffer -----------------
619 * header
620  m_b(ist +1) = 9999999999999d0
621  m_b(ist +2) = 9d12 + id*10 +9d0
622  m_b(ist +3) = 9d12 + iflag1*10 +9d0
623  m_b(ist +4) = 9d12 + iflag2*10 +9d0
624 * dummy vertical scale
625  m_b(ist +5) = -100d0
626  m_b(ist +6) = 100d0
627 * pointer used to speed up search of histogram address
628  m_b(ist +7) = 0d0
629 * information on binning
630  ist2 = ist+7
631  m_b(ist2 +1) = nchx
632  m_b(ist2 +2) = xl
633  m_b(ist2 +3) = xu
634  ddx = xu-xl
635  IF(ddx .EQ. 0d0) CALL glk_stop1('+++GLK_Book1: xl=xu, id= ',id)
636  m_b(ist2 +4) = dfloat(nchx)/ddx
637 *
638 * under/over-flow etc.
639  ist3 = ist+11
640  DO j=1,13
641  m_b(ist3 +j)=0d0
642  ENDDO
643  RETURN
644 *----------------
645  900 CALL glk_retu1(' WARNING GLK_Book1: already exists id= ', id)
646  END
647 
648  SUBROUTINE glk_retu1(mesage,id)
649 * *******************************
650  IMPLICIT NONE
651  include 'GLK.h'
652  SAVE
653  INTEGER id
654  CHARACTER*(*) mesage
655 *-----------------------------
656  WRITE(m_out,'(a)')
657  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
658  WRITE(m_out,'(a,a,i10,a)')
659  $ '++ ', mesage, id, ' ++'
660  WRITE(m_out,'(a)')
661  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
662  WRITE(6 ,'(a)')
663  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
664  WRITE(6 ,'(a,a,i10,a)')
665  $ '++ ', mesage, id, ' ++'
666  WRITE(6 ,'(a)')
667  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
668  END
669 
670  SUBROUTINE glk_stop1(mesage,id)
671 * *******************************
672  IMPLICIT NONE
673  include 'GLK.h'
674  SAVE
675  CHARACTER*(*) mesage
676  INTEGER id
677 *-----------------------------
678  WRITE(m_out,'(a)')
679  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
680  WRITE(m_out,'(a,a,i10,a)')
681  $ '++ ', mesage, id, ' ++'
682  WRITE(m_out,'(a)')
683  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
684  WRITE(6 ,'(a)')
685  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
686  WRITE(6 ,'(a,a,i10,a)')
687  $ '++ ', mesage, id, ' ++'
688  WRITE(6 ,'(a)')
689  $ '++++++++++++++++++++++++++++++++++++++++++++++++++++++'
690  stop
691  END
692 
693 
694  SUBROUTINE glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
695 * ********************************************************
696 * decoding option flags
697 * **********************************
698  IMPLICIT NONE
699  include 'GLK.h'
700  INTEGER id,ioplog,iopsla,ioperb,iopsc1,iopsc2
701  INTEGER ist,iflag1,lact
702 *----------------------------------------------------------------
703  CALL glk_hadres(id,lact)
704  IF(lact .EQ. 0) RETURN
705  ist=m_index(lact,2)
706 * decoding flags
707  iflag1 = nint(m_b(ist+3)-9d0-9d12)/10
708  ioplog = mod(iflag1,10)
709  iopsla = mod(iflag1,100)/10
710  ioperb = mod(iflag1,1000)/100
711  iopsc1 = mod(iflag1,10000)/1000
712  iopsc2 = mod(iflag1,100000)/10000
713  END
714 
715  SUBROUTINE glk_idopt(id,ch)
716 * ************************
717  IMPLICIT NONE
718  include 'GLK.h'
719  INTEGER id
720  CHARACTER*4 ch
721 *
722  INTEGER lact,ist,ioplog,ioperb,iopsla,iopsc1,iopsc2,iflag1
723 *----------------------------------------------------------------
724  CALL glk_hadres(id,lact)
725  IF(lact .EQ. 0) RETURN
726  ist=m_index(lact,2)
727 * decoding flags
728  CALL glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
729  IF(ch .EQ. 'LOGY' ) THEN
730 * log scale for print
731  ioplog = 2
732  ELSEIF(ch .EQ. 'ERRO' ) THEN
733 * errors in printing/plotting
734  ioperb = 2
735  ELSEIF(ch .EQ. 'SLAN' ) THEN
736 * slanted line in plotting
737  iopsla = 2
738  ELSEIF(ch .EQ. 'YMIN' ) THEN
739  iopsc1 = 2
740  ELSEIF(ch .EQ. 'YMAX' ) THEN
741  iopsc2 = 2
742  ENDIF
743 * encoding back
744  iflag1 = ioplog+10*iopsla+100*ioperb+1000*iopsc1+10000*iopsc2
745  m_b(ist+3) = 9d12 + iflag1*10 +9d0
746  END
747 
748 
749  SUBROUTINE glk_bookfun1(id,title,nchx,xmin,xmax,func)
750 */////////////////////////////////////////////////////////////////////////
751 *// fills histogram with function func(x) //
752 */////////////////////////////////////////////////////////////////////////
753  IMPLICIT NONE
754  include 'GLK.h'
755  INTEGER id
756  DOUBLE PRECISION xmin,xmax,func
757  CHARACTER*80 title
758 *
759  DOUBLE PRECISION yy(m_maxnb)
760  EXTERNAL func
761  LOGICAL glk_exist
762  INTEGER ib,nchx
763  DOUBLE PRECISION xl,xu,x
764 *---------------------------------------------------------------------
765  CALL glk_initialize
766  IF(glk_exist(id)) goto 900
767  15 xl=xmin
768  xu=xmax
769  CALL glk_book1(id,title,nchx,xl,xu)
770 *...slanted line in plotting
771  CALL glk_idopt(id,'SLAN')
772  IF(nchx .GT. 200) goto 901
773  DO ib=1,nchx
774  x= xmin +(xmax-xmin)/nchx*(ib-0.5d0)
775  yy(ib) = func(x)
776  ENDDO
777  CALL glk_pak(id,yy)
778  RETURN
779  900 CALL glk_retu1('+++GLK_BookFun1: already exists id=',id)
780  CALL glk_delet(id)
781  goto 15
782  901 CALL glk_stop1('+++GLK_BookFun1: to many bins, id=',id)
783  END
784 
785  SUBROUTINE glk_bookfun1i(id,title,nchx,xmin,xmax,func)
786 */////////////////////////////////////////////////////////////////////////
787 *// Fills histogram with function func(x) //
788 *// Gauss integration over each bin is done, can be slow. //
789 */////////////////////////////////////////////////////////////////////////
790  IMPLICIT NONE
791  include 'GLK.h'
792  INTEGER id
793  DOUBLE PRECISION xmin,xmax,func
794  CHARACTER*80 title
795 *
796  DOUBLE PRECISION yy(m_maxnb)
797  EXTERNAL func
798  LOGICAL glk_exist
799  INTEGER ib,nchx
800  DOUBLE PRECISION xl,xu,x
801  DOUBLE PRECISION glk_gauss,a,b,eeps,dx
802 *---------------------------------------------------------------------
803  CALL glk_initialize
804  IF(glk_exist(id)) goto 900
805  15 xl=xmin
806  xu=xmax
807  CALL glk_book1(id,title,nchx,xl,xu)
808  IF(nchx .GT. 200) goto 901
809  eeps = -0.01d0 !!! relat. precision requirement not very demanding
810  dx = (xmax-xmin)/nchx
811  DO ib=1,nchx
812  a= xmin +dx*(ib-1)
813  b= xmin +dx*ib
814  yy(ib) = glk_gauss(func,a,b,eeps)/dx !! 16-point Gauss integration over bin
815  ENDDO
816  CALL glk_pak(id,yy)
817  RETURN
818  900 CALL glk_retu1('+++GLK_BookFun1I: already exists id=',id)
819  CALL glk_delet(id)
820  goto 15
821  901 CALL glk_stop1('+++GLK_BookFun1I: to many bins, id=',id)
822  END
823 
824  SUBROUTINE glk_bookfun1s(id,title,nchx,xmin,xmax,func)
825 */////////////////////////////////////////////////////////////////////////
826 *// Fills histogram with function func(x) //
827 *// three point fit used //
828 */////////////////////////////////////////////////////////////////////////
829  IMPLICIT NONE
830  include 'GLK.h'
831  DOUBLE PRECISION xmin,xmax,func
832  EXTERNAL func
833  INTEGER id,nchx
834  CHARACTER*80 title
835 * locals
836  DOUBLE PRECISION yy(m_maxnb),yy1(0:m_maxnb)
837  LOGICAL glk_exist
838  DOUBLE PRECISION xl,xu,x3,x2,dx
839  INTEGER ib
840 *--------------------------------------------------------
841  CALL glk_initialize
842  IF( glk_exist(id) ) goto 900
843  15 xl=xmin
844  xu=xmax
845  CALL glk_book1(id,title,nchx,xl,xu)
846 
847 *...slanted line in plotting
848  CALL glk_idopt(id,'SLAN')
849  IF(nchx.gt.200) goto 901
850 
851  yy1(0) = func(xmin)
852  dx=(xmax-xmin)/nchx
853 
854  DO ib=1,nchx
855  x2= xmin +dx*(ib-0.5d0)
856  x3= x2 +dx*0.5d0
857  yy(ib) = func(x2)
858  yy1(ib) = func(x3)
859 *.. simpson
860  yy(ib) = ( yy1(ib-1) +4*yy(ib) +yy1(ib))/6d0
861  ENDDO
862 
863  CALL glk_pak(id,yy)
864  RETURN
865  900 CALL glk_retu1('+++GLK_BookFun1S: already exists, id=',id)
866  CALL glk_delet(id)
867  goto 15
868  901 CALL glk_stop1(' +++GLK_BookFun1S: to many bins, id=',id)
869  END
870 
871  DOUBLE PRECISION FUNCTION glk_gauss(f,a,b,Eeps)
872 *//////////////////////////////////////////////////////////////////////////////
873 *// //
874 *// This is iterative integration procedure //
875 *// originates probably from CERN library //
876 *// it subdivides inegration range until required precision is reached //
877 *// precision is a difference from 8 and 16 point Gauss itegr. result //
878 *// Eeps POSITIVE treated as absolute precision //
879 *// Eeps NEGATIVE treated as relative precision //
880 *// //
881 *//////////////////////////////////////////////////////////////////////////////
882  IMPLICIT NONE
883  DOUBLE PRECISION f,a,b,eeps
884 *
885  DOUBLE PRECISION c1,c2,bb,s8,s16,y,aa,const,delta,eps,u
886  INTEGER i
887 *
888  DOUBLE PRECISION w(12),x(12)
889  EXTERNAL f
890  DATA const /1.0d-19/
891  DATA w
892  1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
893  2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
894  3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
895  4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
896  DATA x
897  1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
898  2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
899  3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
900  4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
901 *-----------------------------------------------------------------------------
902  eps=abs(eeps)
903  delta=const*abs(a-b)
904  glk_gauss=0d0
905  aa=a
906  5 y=b-aa
907  IF(abs(y) .LE. delta) RETURN
908  2 bb=aa+y
909  c1=0.5d0*(aa+bb)
910  c2=c1-aa
911  s8=0d0
912  s16=0d0
913  DO 1 i=1,4
914  u=x(i)*c2
915  1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
916  DO 3 i=5,12
917  u=x(i)*c2
918  3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
919  s8=s8*c2
920  s16=s16*c2
921  IF(eeps .LT. 0d0) THEN
922  IF(abs(s16-s8) .GT. eps*abs(s16)) goto 4
923  ELSE
924  IF(abs(s16-s8) .GT. eps) goto 4
925  ENDIF
926  glk_gauss=glk_gauss+s16
927  aa=bb
928  goto 5
929  4 y=0.5d0*y
930  IF(abs(y) .GT. delta) goto 2
931  WRITE(*,7)
932  glk_gauss=0d0
933  RETURN
934  7 FORMAT(1x,36hgaus ... too high accuracy required)
935  END
936 
937 
938 
939  SUBROUTINE glk_book2(ID,TITLE,NCHX,XL,XU,NCHY,YL,YU)
940 * ***************************************************
941  IMPLICIT NONE
942  include 'GLK.h'
943  INTEGER id,nchx,nchy
944  DOUBLE PRECISION xl,xu,yl,yu
945  CHARACTER*80 title
946 *
947  INTEGER ist,lact,lengt2,j,nnchx,nnchy
948  LOGICAL glk_exist
949 *-------------------------------------------------------------------------
950  CALL glk_initialize
951  IF(glk_exist(id)) goto 900
952  ist=m_length
953  CALL glk_hadres(0,lact)
954  IF(lact .EQ. 0) goto 901
955  m_index(lact,1)=id
956  m_index(lact,2)=m_length
957  CALL glk_copch(title,m_titlc(lact))
958  nnchx=nchx
959  nnchy=nchy
960  lengt2 = m_length +44+nnchx*nnchy
961  IF(lengt2 .GE. m_lenmb) goto 902
962  DO 10 j=m_length+1,lengt2+1
963  10 m_b(j) = 0d0
964  m_length=lengt2
965  m_b(ist+1)=nnchx
966  m_b(ist+2)=xl
967  m_b(ist+3)=xu
968  m_b(ist+4)=float(nnchx)/(m_b(ist+3)-m_b(ist+2))
969  m_b(ist+5)=nnchy
970  m_b(ist+6)=yl
971  m_b(ist+7)=yu
972  m_b(ist+8)=float(nnchy)/(m_b(ist+7)-m_b(ist+6))
973  RETURN
974  900 CALL glk_retu1('GLK_Book2: histo already exists!!!! id=',id)
975  RETURN
976  901 CALL glk_stop1('GLK_Book2: too many histos !!!!! lact= ',lact)
977  RETURN
978  902 CALL glk_stop1('GLK_Book2: too litle storage, m_LenmB=',m_lenmb)
979  RETURN
980  END
981 
982  SUBROUTINE glk_printall
983 * ***********************
984  IMPLICIT NONE
985  include 'GLK.h'
986  SAVE
987  INTEGER i,id
988 
989  DO i=1,m_idmax
990  id=m_index(i,1)
991  IF(id .GT. 0) CALL glk_print(id)
992  ENDDO
993  END
994 
995  SUBROUTINE glk_listprint(mout)
996 *//////////////////////////////////////////////////////////////////////////////
997 *// //
998 *//////////////////////////////////////////////////////////////////////////////
999  IMPLICIT NONE
1000  include 'GLK.h'
1001  INTEGER i,id
1002  CHARACTER*80 title
1003  INTEGER nb,mout
1004  DOUBLE PRECISION xmin,xmax
1005 *----------------------------------
1006  WRITE(mout,*)
1007  $'============================================================================================'
1008  WRITE(mout,*)
1009  $' ID TITLE nb, xmin, xmax'
1010  DO i=1,m_idmax
1011  id=m_index(i,1)
1012  IF(id .NE. 0) THEN
1013  CALL glk_hinbo1(id,title,nb,xmin,xmax)
1014  WRITE(mout,'(i8,a,a,i8,2g14.6)') id, ' ', title, nb,xmin,xmax
1015  ENDIF
1016  ENDDO
1017  END
1018 
1019 
1020 
1021  SUBROUTINE glk_print(id)
1022 * ***********************
1023  IMPLICIT NONE
1024  include 'GLK.h'
1025  INTEGER id
1026 *
1027  DOUBLE PRECISION xl,bind,xlow,z,er,avex,dx,fact,ovef,undf,bmax,bmin,deltb
1028  DOUBLE PRECISION sum,sumw,sumx
1029  INTEGER ist,ist2,ist3,idec,k2,k1,kros,j,ind,i,n,i1,ky,nchy,kx,nent,iflag2,lmx
1030  INTEGER ioplog,iopsla,ioperb,iopsc1,iopsc2,lact,ker,ityphi,kzer,k,ibn,nchx,istr
1031  LOGICAL llg
1032  CHARACTER*1 line(0:105),lchr(22),lb,lx,li,l0
1033  SAVE lb,lx,li,l0,lchr
1034  DATA lb,lx,li,l0 /' ','X','I','0'/
1035  DATA lchr/' ','1','2','3','4','5','6','7','8','9',
1036  $ 'A','B','C','D','E','F','G','H','I','J','K','*'/
1037 *---------------------------------------------------------------------------------
1038  CALL glk_hadres(id,lact)
1039  IF(lact .EQ. 0) goto 900
1040  ist = m_index(lact,2)
1041  ist2 = ist+7
1042  ist3 = ist+11
1043  idec = nint(m_b(ist+2)-9d0-9d12)/10
1044  IF(idec .NE. id) WRITE(6,*) '++++GLK_PRINT: PANIC! ID,IDEC= ',id,idec
1045  CALL glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
1046  ker = ioperb-1
1047  lmx = 67
1048  IF(ker .EQ. 1) lmx=54
1049  nent=m_index(lact,3)
1050  IF(nent .EQ. 0) goto 901
1051  WRITE(m_out,1000) id,m_titlc(lact)
1052  1000 FORMAT('1',/,1x,i9,10x,a)
1053 *
1054 * one-dim. histo
1055  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1056  ityphi = mod(iflag2,10)
1057  IF(ityphi .EQ. 2) goto 200
1058  IF( (ityphi.NE.1) .AND. (ityphi.NE.3) )
1059  $ CALL glk_stop1(' GLK_PRINT: wrong histo type, id=',id)
1060 
1061  nchx = m_b(ist2 +1)
1062  xl = m_b(ist2 +2)
1063  dx = ( m_b(ist2 +3)-m_b(ist2 +2) )/float(nchx)
1064 * fixing vertical scale
1065  istr=ist+m_buf1+1
1066  bmin = m_b(istr)
1067  bmax = m_b(istr)+1d-5*abs(m_b(istr)) ! problems for single bin case
1068  DO ibn=istr,istr+nchx-1
1069  bmax = max(bmax,m_b(ibn))
1070  bmin = min(bmin,m_b(ibn))
1071  ENDDO
1072  IF(bmin .EQ. bmax) goto 903
1073  IF(iopsc1 .EQ. 2) bmin=m_b(ist +5)
1074  IF(iopsc2 .EQ. 2) bmax=m_b(ist +6)
1075 *
1076  llg=ioplog .EQ. 2
1077  IF(llg.and.bmin .LE. 0d0) bmin=bmax/10000.d0
1078 *
1079  deltb = bmax-bmin
1080  IF(deltb .EQ. 0d0) goto 902
1081  fact = (lmx-1)/deltb
1082  kzer = -bmin*fact+1.00001d0
1083  IF(llg) fact=(lmx-1)/(log(bmax)-log(bmin))
1084  IF(llg) kzer=-log(bmin)*fact+1.00001d0
1085 *
1086  undf = m_b(ist3 +1)
1087  ovef = m_b(ist3 +3)
1088  avex = 0d0
1089  sumw = m_b(ist3 +8)
1090  sumx = m_b(ist3 +9)
1091  IF(sumw .NE. 0d0) avex = sumx/sumw
1092  WRITE(m_out,'(4a15 )') 'nent',' sum','bmin','bmax'
1093  WRITE(m_out,'(i15,3e15.5)') nent, sum, bmin, bmax
1094  WRITE(m_out,'(4a15 )') 'undf','ovef','sumw','avex'
1095  WRITE(m_out,'(4e15.5)') undf, ovef, sumw, avex
1096 *
1097  IF(llg) write(m_out,1105)
1098  1105 format(35x,17hlogarithmic scale)
1099 *
1100  kzer=max0(kzer,0)
1101  kzer=min0(kzer,lmx)
1102  xlow=xl
1103  do 100 k=1,nchx
1104 * first fill with blanks
1105  do 45 j=1,105
1106  45 line(j) =lb
1107 * THEN fill upper and lower boundry
1108  line(1) =li
1109  line(lmx)=li
1110  ind=istr+k-1
1111  bind=m_b(ind)
1112  bind= max(bind,bmin)
1113  bind= min(bind,bmax)
1114  kros=(bind-bmin)*fact+1.0001d0
1115  IF(llg) kros=log(bind/bmin)*fact+1.0001d0
1116  k2=max0(kros,kzer)
1117  k2=min0(lmx,max0(1,k2))
1118  k1=min0(kros,kzer)
1119  k1=min0(lmx,max0(1,k1))
1120  do 50 j=k1,k2
1121  50 line(j)=lx
1122  line(kzer)=l0
1123  z=m_b(ind)
1124  IF(ker .NE. 1) THEN
1125  WRITE(m_out,'(a, f7.4, a, d14.6, 132a1)')
1126  $ ' ', xlow,' ', z,' ',(line(i),i=1,lmx)
1127  ELSE
1128  er=dsqrt(dabs(m_b(ind+nchx)))
1129  WRITE(m_out,'(a,f7.4, a,d14.6, a,d14.6, 132a1 )')
1130  $ ' ',xlow,' ', z,' ', er,' ',(line(i),i=1,lmx)
1131  ENDIF
1132  xlow=xlow+dx
1133  100 continue
1134  RETURN
1135 *//////////////////////////////////////////////////////////////////////
1136 *// two dimensional requires complete restoration and tests //
1137 *//////////////////////////////////////////////////////////////////////
1138  200 continue
1139  nchx=m_b(ist+1)
1140  nchy=m_b(ist+5)
1141  WRITE(m_out,2000) (lx,i=1,nchy)
1142  2000 format(1h ,10x,2hxx,100a1)
1143  do 300 kx=1,nchx
1144  do 250 ky=1,nchy
1145  k=ist +m_buf2 +kx+nchx*(ky-1)
1146  n=m_b(k)+1.99999d0
1147  n=max0(n,1)
1148  n=min0(n,22)
1149  IF(dabs(m_b(k)) .LT. 1d-20) n=1
1150  line(ky)=lchr(n)
1151  250 continue
1152  line(nchy+1)=lx
1153  i1=nchy+1
1154  WRITE(m_out,2100) (line(i),i=1,i1)
1155  2100 format(1h ,10x,1hx,100a1)
1156  300 continue
1157  WRITE(m_out,2000) (lx,i=1,nchy)
1158  RETURN
1159  900 CALL glk_retu1('GLK_PRINT: nonexisting histo',id)
1160  RETURN
1161  901 CALL glk_retu1(.eq.' GLK_PRINT: nent0',id)
1162  RETURN
1163  902 CALL glk_retu1(' GLK_PRINT: wrong plotting limits, id=',id)
1164  RETURN
1165  903 CALL glk_retu1(.eq.' GLK_PRINT: bminbmax, id=',id)
1166  END
1167 
1168  SUBROUTINE glk_operat(ida,chr,idb,idc,coef1,coef2)
1169 * **********************************************
1170  IMPLICIT NONE
1171  include 'GLK.h'
1172  INTEGER ida,idb,idc
1173  DOUBLE PRECISION coef1,coef2
1174  CHARACTER*80 title
1175  CHARACTER*1 chr
1176 *
1177  DOUBLE PRECISION xl,xu
1178  INTEGER ista,ista2,ista3,ncha,iflag2a,ityphia,lactb
1179  INTEGER k,j,nchc,istc2,istc3,i1,j2,j3,j1,i2,i3,istc,istb2,istb3,nchb
1180  INTEGER lacta,id,istb,nchx,iflag2b,ityphib,lactc
1181 *----------------------------------------------------------
1182  CALL glk_hadres(ida,lacta)
1183  IF(lacta .EQ. 0) RETURN
1184  ista = m_index(lacta,2)
1185  ista2 = ista+7
1186  ista3 = ista+11
1187  ncha = m_b(ista2+1)
1188 * check for type
1189  iflag2a = nint(m_b(ista+4)-9d0-9d12)/10
1190  ityphia = mod(iflag2a,10)
1191  IF(ityphia .NE. 1) CALL glk_stop1('GLK_Operat: 1-dim histos only, id=',id)
1192 *------------------
1193  CALL glk_hadres(idb,lactb)
1194  IF(lactb .EQ. 0) RETURN
1195  istb = m_index(lactb,2)
1196  istb2 = istb+7
1197  istb3 = istb+11
1198  nchb = m_b(istb2+1)
1199  IF(nchb .NE. ncha) goto 900
1200 * check for type
1201  iflag2b = nint(m_b(istb+4)-9d0-9d12)/10
1202  ityphib = mod(iflag2b,10)
1203  IF(ityphib .NE. 1) CALL glk_stop1('GLK_Operat: 1-dim histos only, id=',id)
1204 *------------------
1205  CALL glk_hadres(idc,lactc)
1206  IF(lactc .EQ. 0) THEN
1207 * ...if nonexistent, histo idc is here defined
1208  CALL glk_hinbo1(ida,title,nchx,xl,xu)
1209  CALL glk_book1(idc,title,nchx,xl,xu)
1210  CALL glk_hadres(idc,lactc)
1211  istc = m_index(lactc,2)
1212 *...option copied from ida
1213  m_b(istc+ 3)= m_b(ista +3)
1214  ENDIF
1215 *...one nominal entry recorded
1216  m_index(lactc,3) = 1
1217 *
1218  istc = m_index(lactc,2)
1219  istc2 = istc+7
1220  istc3 = istc+11
1221  nchc = m_b(istc2+1)
1222 *
1223  IF(nchc .NE. ncha) goto 900
1224  IF(ncha .NE. nchb .OR. nchb .NE. nchc) goto 900
1225  DO k=1,ncha+3
1226  IF(k .GT. ncha) THEN
1227 * underflow, overflow
1228  j=k-ncha
1229  i1 = ista3 +j
1230  i2 = istb3 +j
1231  i3 = istc3 +j
1232  j1 = ista3 +3+j
1233  j2 = istb3 +3+j
1234  j3 = istc3 +3+j
1235  ELSE
1236 * normal bins
1237  i1 = ista +m_buf1 +k
1238  i2 = istb +m_buf1 +k
1239  i3 = istc +m_buf1 +k
1240  j1 = ista +m_buf1 +ncha+k
1241  j2 = istb +m_buf1 +ncha+k
1242  j3 = istc +m_buf1 +ncha+k
1243  ENDIF
1244  IF (chr .EQ. '+') THEN
1245  m_b(i3) = coef1*m_b(i1) + coef2*m_b(i2)
1246  m_b(j3) = coef1**2*m_b(j1) + coef2**2*m_b(j2)
1247  ELSEIF(chr .EQ. '-') THEN
1248  m_b(i3) = coef1*m_b(i1) - coef2*m_b(i2)
1249  m_b(j3) = coef1**2*m_b(j1) + coef2**2*m_b(j2)
1250  ELSEIF(chr .EQ. '*') THEN
1251  m_b(j3) = (coef1*coef2)**2
1252  $ *(m_b(j1)*m_b(i2)**2 + m_b(j2)*m_b(i1)**2)
1253  m_b(i3) = coef1*m_b(i1) * coef2*m_b(i2)
1254  ELSEIF(chr .EQ. '/') THEN
1255  IF(m_b(i2) .EQ. 0d0) THEN
1256  m_b(i3) = 0d0
1257  m_b(j3) = 0d0
1258  ELSE
1259 *** m_b(j3) = (coef1/coef2)**2/m_b(i2)**4 ! problems with overflow
1260 *** $ *(m_b(j1)*m_b(i2)**2 + m_b(j2)*m_b(i1)**2) ! problems with overflow
1261  m_b(j3) = (coef1/coef2)**2 *m_b(j1) /m_b(i2)**2
1262  $ +(coef1/coef2)**2 *m_b(j2) *(m_b(i1)/m_b(i2)**2)**2
1263  m_b(i3) = (coef1*m_b(i1) )/( coef2*m_b(i2))
1264  ENDIF
1265  ELSE
1266  goto 901
1267  ENDIF
1268  ENDDO
1269  RETURN
1270  900 WRITE(m_out,*) '+++++ GLK_Operat: non-equal no. bins ',ida,idb,idc
1271  WRITE( 6,*) '+++++ GLK_Operat: non-equal no. bins ',ida,idb,idc
1272  stop
1273  901 WRITE(m_out,*) '+++++ GLK_Operat: wrong chr=',chr
1274  WRITE( 6,*) '+++++ GLK_Operat: wrong chr=',chr
1275  stop
1276  END
1277 
1278  SUBROUTINE glk_hinbo1(id,title,nchx,xl,xu)
1279 * **************************************
1280  IMPLICIT NONE
1281  include 'GLK.h'
1282  INTEGER id,nchx
1283  DOUBLE PRECISION xl,xu
1284  CHARACTER*80 title
1285  INTEGER lact,ist,ist2
1286 *----------------------------------------------------------------------
1287  CALL glk_hadres(id,lact)
1288  IF(lact .EQ. 0) THEN
1289  CALL glk_stop1('+++STOP in GLK_hinbo1: wrong id=',id)
1290  ENDIF
1291  ist=m_index(lact,2)
1292  ist2 = ist+7
1293  nchx = m_b(ist2 +1)
1294  xl = m_b(ist2 +2)
1295  xu = m_b(ist2 +3)
1296  title = m_titlc(lact)
1297  END
1298 
1299  SUBROUTINE glk_unpak(id,a,chd1,idum)
1300 * *********************************
1301 * getting out histogram content (and error)
1302 * chd1= 'ERRO' is nonstandard option (unpack errors)
1303 * ***********************************
1304  IMPLICIT NONE
1305  include 'GLK.h'
1306  INTEGER id,idum
1307  DOUBLE PRECISION a(*)
1308  CHARACTER*(*) chd1
1309 *
1310  INTEGER lact,ist,ist2,iflag2,ityphi,local,nch,nchy,ib
1311 *------------------------------------------------------------------------
1312  CALL glk_hadres(id,lact)
1313  IF(lact .EQ. 0) goto 900
1314  ist = m_index(lact,2)
1315  ist2 = ist+7
1316  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1317  ityphi = mod(iflag2,10)
1318  IF(ityphi .EQ. 1) THEN
1319  nch = m_b(ist2 +1)
1320  local = ist +m_buf1
1321  ELSEIF(ityphi .EQ. 2) THEN
1322  nchy = m_b(ist2+5)
1323  nch = nch*nchy
1324  local = ist+ m_buf2
1325  ELSE
1326  CALL glk_stop1('+++GLK_UnPak: check type of histo id=',id)
1327  ENDIF
1328  do 10 ib=1,nch
1329  IF(chd1 .NE. 'ERRO') THEN
1330 * normal bin
1331  a(ib) = m_b(local+ib)
1332  ELSE
1333 * error content
1334  IF(ityphi .EQ. 2) goto 901
1335  a(ib) = dsqrt( dabs(m_b(local+nch+ib) ))
1336  ENDIF
1337  10 continue
1338  RETURN
1339  900 CALL glk_retu1('+++GLK_UnPak: nonexisting id=',id)
1340  RETURN
1341  901 CALL glk_retu1('+++GLK_UnPak: no errors, two-dim, id=',id)
1342  END
1343 
1344  SUBROUTINE glk_pak(id,a)
1345 * ************************
1346 * Loading in histogram content
1347 * ***********************************
1348  IMPLICIT NONE
1349  include 'GLK.h'
1350  INTEGER id
1351  DOUBLE PRECISION a(*)
1352 *
1353  INTEGER lact,ist,ist2,iflag2,ityphi,nch,local,ib,nchy
1354 *----------------------------------------------------
1355  CALL glk_hadres(id,lact)
1356  IF(lact .EQ. 0) goto 900
1357  ist = m_index(lact,2)
1358  ist2 = ist+7
1359 * 2-dimens histo alowed
1360  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1361  ityphi = mod(iflag2,10)
1362  IF(ityphi .EQ. 1) THEN
1363  nch = m_b(ist2 +1)
1364  local = ist+m_buf1
1365  ELSEIF(ityphi .EQ. 2) THEN
1366  nchy = m_b(ist2+5)
1367  nch = nch*nchy
1368  local = ist+m_buf2
1369  ELSE
1370  CALL glk_stop1('+++GLK_Pak: wrong histo type, id=',id)
1371  ENDIF
1372  DO ib=1,nch
1373  m_b(local +ib) = a(ib)
1374  ENDDO
1375 * one nominal entry recorded
1376  m_index(lact,3) = 1
1377  RETURN
1378  900 CONTINUE
1379  CALL glk_stop1('+++GLK_Pak: nonexisting id=',id)
1380  END
1381 
1382  SUBROUTINE glk_pake(id,a)
1383 * **********************
1384 * Loading in error content
1385 * ***********************************
1386  IMPLICIT NONE
1387  include 'GLK.h'
1388  INTEGER id
1389  DOUBLE PRECISION a(*)
1390 *
1391  INTEGER lact,ist,ist2,nch,iflag2,ityphi
1392  INTEGER nb,ib
1393 *---------------------------------------------------------------------
1394  CALL glk_hadres(id,lact)
1395  IF(lact .EQ. 0) goto 901
1396  ist = m_index(lact,2)
1397  ist2 = ist+7
1398  nch=m_b(ist2+1)
1399 * 2-dimens histo NOT alowed
1400  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1401  ityphi = mod(iflag2,10)
1402  IF(ityphi .NE. 1) goto 900
1403  DO ib=1,nch
1404  m_b(ist+m_buf1+nch+ib) = a(ib)**2
1405  ENDDO
1406  CALL glk_idopt( id,'ERRO')
1407  RETURN
1408  900 CALL glk_stop1('+++GLK_Pake: only for 1-dim hist, id=',id)
1409  RETURN
1410  901 CALL glk_stop1('+++GLK_Pake: nonexisting id=',id)
1411  END
1412 
1413 
1414  SUBROUTINE glk_range1(id,ylr,yur)
1415 * *****************************
1416 * provides y-scale for 1-dim plots
1417 * ***********************************
1418  IMPLICIT NONE
1419  include 'GLK.h'
1420  INTEGER id
1421  DOUBLE PRECISION ylr,yur
1422 *
1423  INTEGER lact,ist,ist2,nch,ib,ioplog,iopsla,ioperb,iopsc1,iopsc2
1424  DOUBLE PRECISION yl,yu
1425 *-------------------------------------------------------------
1426  CALL glk_hadres(id,lact)
1427  IF(lact .EQ. 0) RETURN
1428  ist = m_index(lact,2)
1429  ist2 = ist+7
1430  nch = m_b(ist2 +1)
1431  yl = m_b(ist+m_buf1+1)
1432  yu = m_b(ist+m_buf1+1)
1433  DO ib=1,nch
1434  yl = min(yl,m_b(ist+m_buf1+ib))
1435  yu = max(yu,m_b(ist+m_buf1+ib))
1436  ENDDO
1437 * For default range some safety range is added
1438  yu = yu + 0.05*abs(yu-yl)
1439 *** yl = yl - 0.05*ABS(yu-yl) ! to be decided later on
1440 
1441 * If range was pre-defined then yl,yu are overwritten
1442  CALL glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
1443  IF(iopsc1 .EQ. 2) yl= m_b( ist +5)
1444  IF(iopsc2 .EQ. 2) yu= m_b( ist +6)
1445  ylr = yl
1446  yur = yu
1447  END
1448 
1449 
1450  SUBROUTINE glk_hinbo2(id,nchx,xl,xu,nchy,yl,yu)
1451 * *******************************************
1452  IMPLICIT NONE
1453  include 'GLK.h'
1454  INTEGER id,nchx,nchy
1455  DOUBLE PRECISION xl,xu,yl,yu
1456  INTEGER lact,ist,ist2
1457 *--------------------------------------------------
1458  CALL glk_hadres(id,lact)
1459  IF(lact .EQ. 0) goto 900
1460  ist = m_index(lact,2)
1461  ist2 = ist+7
1462  nchx = m_b(ist2 +1)
1463  xl = m_b(ist2 +2)
1464  xu = m_b(ist2 +3)
1465  nchy = m_b(ist2 +5)
1466  yl = m_b(ist2 +6)
1467  yu = m_b(ist2 +7)
1468  RETURN
1469  900 CALL glk_stop1(' +++GLK_hinbo2: nonexisting histo id= ',id)
1470  END
1471 
1472 
1473  SUBROUTINE glk_ymaxim(id,wmax)
1474 * **************************
1475  IMPLICIT NONE
1476  include 'GLK.h'
1477  INTEGER id
1478  DOUBLE PRECISION wmax
1479  INTEGER lact,ist,jd,k
1480 *-------------------------------------------------------
1481  IF(id .NE. 0) THEN
1482  CALL glk_hadres(id,lact)
1483  IF(lact .EQ. 0) RETURN
1484  ist= m_index(lact,2)
1485  m_b(ist+6) =wmax
1486  CALL glk_idopt(id,'YMAX')
1487  ELSE
1488  DO k=1,m_idmax
1489  IF(m_index(k,1) .EQ. 0) goto 20
1490  ist=m_index(k,2)
1491  jd =m_index(k,1)
1492  m_b(ist+6) =wmax
1493  CALL glk_idopt(jd,'YMAX')
1494  ENDDO
1495  20 CONTINUE
1496  ENDIF
1497  END
1498 
1499  SUBROUTINE glk_yminim(id,wmin)
1500 * ******************************
1501  IMPLICIT NONE
1502  include 'GLK.h'
1503  INTEGER id
1504  DOUBLE PRECISION wmin
1505  INTEGER lact,ist,k,jd
1506 *---------------------------------------------
1507  IF(id .NE. 0) THEN
1508  CALL glk_hadres(id,lact)
1509  IF(lact .EQ. 0) RETURN
1510  ist =m_index(lact,2)
1511  m_b(ist+5) =wmin
1512  CALL glk_idopt(id,'YMIN')
1513  ELSE
1514  DO k=1,m_idmax
1515  IF(m_index(k,1) .EQ. 0) goto 20
1516  ist=m_index(k,2)
1517  jd =m_index(k,1)
1518  m_b(ist+5) =wmin
1519  CALL glk_idopt(jd,'YMIN')
1520  ENDDO
1521  20 CONTINUE
1522  ENDIF
1523  END
1524 
1525  SUBROUTINE glk_reset(id,chd1)
1526 * **************************
1527  IMPLICIT NONE
1528  include 'GLK.h'
1529  INTEGER id
1530  CHARACTER*(*) chd1
1531  INTEGER lact,ist,ist2,iflag2,ityphi,ist3,nchx,nch,local,nchy,j
1532 *-------------------------------------------
1533  CALL glk_hadres(id,lact)
1534  IF(lact .LE. 0) RETURN
1535  ist =m_index(lact,2)
1536  ist2 = ist+7
1537 *
1538  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1539  ityphi = mod(iflag2,10)
1540  IF(ityphi .EQ. 1) THEN
1541 * one-dim.
1542  ist3 = ist+11
1543  nchx = m_b(ist2 +1)
1544  nch = 2*nchx
1545  local = ist + m_buf1
1546  ELSEIF(ityphi .EQ. 2) THEN
1547 * two-dim.
1548  ist3 = ist+15
1549  nchx = m_b(ist2 +1)
1550  nchy = m_b(ist2 +5)
1551  nch = nchx*nchy
1552  local = ist +m_buf2
1553  ELSE
1554  CALL glk_stop1('+++GLK_Reset: wrong type id=',id)
1555  ENDIF
1556 * reset miscaelaneous entries and bins
1557  DO j=ist3+1,local +nch
1558  m_b(j) = 0d0
1559  ENDDO
1560 * and no. of entries in m_index
1561  m_index(lact,3) = 0
1562  END
1563 
1564  SUBROUTINE glk_delet(id1)
1565 * *********************
1566 * Now it should work (stj Nov. 91) but watch out!
1567 * should works for 2-dim histos, please check this!
1568 * ***********************************
1569  IMPLICIT NONE
1570  include 'GLK.h'
1571  INTEGER id1
1572 *
1573  LOGICAL glk_exist
1574  INTEGER id,lact,ist,ist2,nch,iflag2,ityphi,local,k,i,l,next,idec,nchx,nchy
1575 *--------------------------------------------
1576  id=id1
1577  IF(id .EQ. 0) goto 300
1578  IF( .NOT. glk_exist(id)) goto 900
1579  CALL glk_hadres(id,lact)
1580  ist = m_index(lact,2)
1581  ist2 = ist+7
1582 *----
1583 *[[[ WRITE(6,*) 'GLK_DELET-ing ID= ',ID
1584  idec = nint(m_b(ist+2)-9d0-9d12)/10
1585  IF(idec .NE. id) WRITE(6,*)
1586  $ '++++GLK_DELET: ALARM! ID,IDEC= ',id,idec
1587 *----
1588  nch = m_b(ist2 +1)
1589  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
1590  ityphi = mod(iflag2,10)
1591  IF(ityphi .EQ. 1) THEN
1592 * one-dim.
1593  nchx = m_b(ist2 +1)
1594  nch = 2*nchx
1595 * lenght of local histo to be removed
1596  local = nch+m_buf1+1
1597  ELSEIF(ityphi .EQ. 2) THEN
1598 * two-dim.
1599  nchx = m_b(ist2 +1)
1600  nchy = m_b(ist2 +5)
1601  nch = nchx*nchy
1602 * lenght of local histo to be removed
1603  local = nch+m_buf2+1
1604  ELSE
1605  CALL glk_stop1('+++GLK_Delet: wrong type id=',id)
1606  ENDIF
1607 * starting position of next histo in storage b
1608  next = ist+1 +local
1609 * move down all histos above this one
1610  DO 15 k =next,m_length
1611  m_b(k-local)=m_b(k)
1612  15 CONTINUE
1613 * define new end of storage
1614  m_length=m_length-local
1615 * clean free space at the end of storage b
1616  DO 20 k=m_length+1, m_length+local
1617  20 m_b(k)=0d0
1618 * shift adresses of all displaced histos
1619  DO 25 l=lact+1,m_idmax
1620  IF(m_index(l,1) .NE. 0) m_index(l,2)=m_index(l,2)-local
1621  25 CONTINUE
1622 * move entries in m_index down by one and remove id=lact entry
1623  DO 30 l=lact+1,m_idmax
1624  m_index(l-1,1)=m_index(l,1)
1625  m_index(l-1,2)=m_index(l,2)
1626  m_index(l-1,3)=m_index(l,3)
1627  m_titlc(l-1)=m_titlc(l)
1628  30 CONTINUE
1629 * last entry should be always empty
1630  m_index(m_idmax,1)=0
1631  m_index(m_idmax,2)=0
1632  m_index(m_idmax,3)=0
1633  do 50 k=1,80
1634  50 m_titlc(m_idmax)(k:k)=' '
1635  RETURN
1636 * -----------------------------------
1637 * Deleting all histos at once!!!
1638  300 m_length=0
1639  DO 400 i=1,m_idmax
1640  DO 340 k=1,3
1641  340 m_index(i,k)=0
1642  DO 350 k=1,80
1643  350 m_titlc(i)(k:k)=' '
1644  400 CONTINUE
1645  RETURN
1646 * -----------------------------------
1647  900 CONTINUE
1648  CALL glk_retu1(' +++GLK_DELET: nonexisting histo id= ',id)
1649  END
1650 
1651 
1652  SUBROUTINE glk_copch(ch1,ch2)
1653 * *****************************
1654  IMPLICIT NONE
1655 * copies CHARACTER*80 ch1 into ch2 up to a first $ sign
1656  CHARACTER*80 ch1,ch2
1657  LOGICAL met
1658  INTEGER i
1659 *----------------------------
1660  met = .false.
1661  DO i=1,80
1662  IF( ch1(i:i) .EQ. '$' .or. met ) THEN
1663  ch2(i:i)=' '
1664  met=.true.
1665  ELSE
1666  ch2(i:i)=ch1(i:i)
1667  ENDIF
1668  ENDDO
1669  END
1670 
1671  INTEGER FUNCTION glk_jadre2(id)
1672 *------------------------------------------------
1673 * Good old version -- but it is very very slow!!!
1674 * In the case of 100 histograms or more.
1675 *------------------------------------------------
1676  IMPLICIT NONE
1677  include 'GLK.h'
1678  INTEGER id,i
1679 *---------------------------------------
1680  glk_jadre2=0
1681  DO 1 i=1,m_idmax
1682  IF(m_index(i,1) .EQ. id) goto 2
1683  1 CONTINUE
1684 * Nothing found.
1685  RETURN
1686 * Found: id=0 is also legitimate find!!!
1687  2 glk_jadre2=i
1688  END
1689 
1690  SUBROUTINE glk_hadres(id1,jadres)
1691 * *********************************
1692 *--------------------------------------------------------------------
1693 * Educated guess based on past history is used to find quickly
1694 * location of the histogram in the matrix m_index.
1695 * This is based on observation that subsequent histogram calls
1696 * are linked into loops (so one can predict easily which histo will
1697 * be called next time).
1698 *--------------------------------------------------------------------
1699  IMPLICIT NONE
1700  include 'GLK.h'
1701  INTEGER id1,jadres
1702  INTEGER ist,i,id
1703 *----------------------------------------------------------------------
1704  INTEGER iguess,jdlast,idlast
1705  DATA iguess,jdlast,idlast /-2141593,-3141593,-3141593/
1706  SAVE iguess,jdlast,idlast
1707 *----------------------------------------------------------------------
1708  id=id1
1709 * --- The case of ID=0 treated separately, it is used to find out
1710 * --- last entry in the m_index (it is marked with zero)
1711  IF(id .EQ. 0) THEN
1712  DO i=1,m_idmax
1713  IF(m_index(i,1) .EQ. 0) goto 44
1714  ENDDO
1715  CALL glk_stop1('+++GLK_hadres: Too short m_index=',m_index)
1716  44 CONTINUE
1717  jadres = i
1718  RETURN
1719  ENDIF
1720 
1721 * --- Omit sophistications if lack of initialization
1722  IF(jdlast .EQ. -3141593) goto 10
1723  IF(iguess .EQ. -2141593) goto 10
1724  IF(iguess .EQ. 0) goto 10
1725  IF(jdlast .EQ. 0) goto 10
1726 
1727 * --- Try first previous histo (for repeated calls)
1728  IF(jdlast .LT. 1 .OR. jdlast .GT. m_idmax) THEN
1729  WRITE(6,*) '+++++ hadres: jdlast=',jdlast
1730  ENDIF
1731  IF(m_index(jdlast,1) .EQ. id) THEN
1732  jadres = jdlast
1733 *## WRITE(6,*)
1734 *## $ 'found, guess based on previous call to jadres ',jdlast
1735  goto 20
1736  ENDIF
1737 
1738 * --- Try current guess based on previous call
1739  IF(iguess .LT. 1 .OR. iguess .GT. m_idmax) THEN
1740  WRITE(6,*)'+++++ hadres: iguess=',iguess
1741  ENDIF
1742  IF(m_index(iguess,1) .EQ. id) THEN
1743  jadres = iguess
1744 *## WRITE(6,*)
1745 *## $ 'found, guess on previous calls recorded in m_b(ist+7)',jdlast
1746  goto 20
1747  ENDIF
1748 
1749 * ================================================
1750 * Do it HARD WAY, Search all matrix m_index
1751 * ================================================
1752  10 CONTINUE
1753 *## WRITE(6,*) 'trying HARD WAY'
1754  DO i=1,m_idmax
1755  jadres=i
1756  IF(m_index(i,1) .EQ. id) goto 20
1757  ENDDO
1758 * -------------------------------------
1759 * Nothing found: jadres=0
1760 * -------------------------------------
1761  jadres=0
1762  RETURN
1763 * =====================================
1764 * Found: Set new guess for next call
1765 * =====================================
1766  20 CONTINUE
1767 * --- and store result as a new guess in previous histo
1768 * --- but only if it existed!!!!
1769  DO i=1,m_idmax
1770  IF(m_index(i,1) .EQ. 0) goto 40
1771  IF(m_index(i,1) .EQ. idlast) THEN
1772  ist=m_index(i,2)
1773  IF(ist .GT. 0 .AND. ist .LT. m_lenmb) m_b(ist +7) = jadres
1774 *## WRITE(6,*) 'STORED id=',id
1775  goto 40
1776  ENDIF
1777  ENDDO
1778  40 CONTINUE
1779 *## WRITE(6,*) 'found, hard way searching all of m_index)', jdlast
1780  iguess = m_b( m_index(jadres,2) +7)
1781  jdlast = jadres
1782  idlast = id
1783  END
1784 
1785 
1786 *--------------------------------------------------------------
1787 * ----------- storing histograms in the disk file -------------
1788 *--------------------------------------------------------------
1789 
1790  SUBROUTINE glk_readfile(Dname)
1791 * ******************************
1792  IMPLICIT NONE
1793  include 'GLK.h'
1794  SAVE
1795  INTEGER ninph
1796  CHARACTER*60 dname
1797 *-------------------------------------------------
1798  CALL glk_initialize
1799 * Read histograms
1800  WRITE( *,*) 'GLK_ReadFile: Reading histos from ', dname
1801  WRITE(m_out,*) 'GLK_ReadFile: Reading histos from ', dname
1802  ninph = 21
1803  OPEN(ninph,file=dname) ! Open dump-file
1804  CALL glk_hrfile(ninph,' ',' ') ! Define dump-file unit in GKL
1805  CALL glk_hrin(0,0,0) ! Read histos from dump-file
1806  CALL glk_hrend(' ') ! Close dump-file
1807  END
1808 
1809  SUBROUTINE glk_writefile(Dname)
1810 * ******************************
1811  IMPLICIT NONE
1812  include 'GLK.h'
1813  SAVE
1814  INTEGER nouth
1815  CHARACTER*60 dname
1816 *-------------------------------------------------
1817  CALL glk_initialize
1818 * Write all histograms into disk file
1819  WRITE( *,*) 'GLK_WriteFile: Writing histos into ', dname
1820  WRITE(m_out,*) 'GLK_WriteFile: Writing histos into ', dname
1821  nouth=22
1822  OPEN(nouth,file=dname) ! Open dump-file
1823  CALL glk_hrfile(nouth,' ',' ') ! Define dump-file in GLK
1824  CALL glk_hrout(0,0,' ') ! Dump all histos on disk
1825  CALL glk_hrend(' ') ! Close dump-file
1826  END
1827 
1828  SUBROUTINE glk_hrfile(nhruni,chd1,chd2)
1829 * ***************************************
1830  IMPLICIT NONE
1831  CHARACTER*(*) chd1,chd2
1832  include 'GLK.h'
1833  SAVE
1834  INTEGER nhruni
1835 *-------------------------------------------------
1836  CALL glk_initialize
1837  m_huni=nhruni
1838  END
1839 
1840  SUBROUTINE glk_hrout(idum1,idum2,chdum)
1841 * ***************************************
1842  IMPLICIT NONE
1843  CHARACTER*8 chdum
1844 *
1845  include 'GLK.h'
1846  SAVE
1847  INTEGER i,k,last,idum1,idum2
1848 *-----------------------------------------------------------------------
1849  CALL glk_initialize
1850  CALL glk_hadres(0,last)
1851  IF(last.EQ.0) goto 900
1852  m_lenind = last -1
1853  WRITE(m_huni,'(6i10)') m_version, m_lenind, m_lenmb, m_length
1854  WRITE(m_huni,'(6i10)') ((m_index(i,k),k=1,3),i=1,m_lenind)
1855  WRITE(m_huni,'(a80)') (m_titlc(i), i=1,m_lenind)
1856  WRITE(m_huni,'(3d24.16)') (m_b(i), i=1,m_length)
1857  RETURN
1858  900 CONTINUE
1859  WRITE(m_out,*) '+++ GLK_hrout: no space in index'
1860  WRITE( *,*) '+++ GLK_hrout: no space in index'
1861  END
1862 
1863 
1864  SUBROUTINE glk_hrin(idum1,idum2,idum3)
1865 * **************************************
1866 * New version which has a possibility to
1867 * MERGE histograms
1868 * If given ID already exists then it is modified by adding 1000000 !!!!
1869 * Mergigng is done simply by appending new histograms at the
1870 * very end of the m_index and bin matrices.
1871 * ***********************************
1872  IMPLICIT NONE
1873  include 'GLK.h'
1874  INTEGER idum1,idum2,idum3
1875  INTEGER l,lact,lenold,istn,idn,k,lenind3,nvrs3,nouth
1876  INTEGER i,lengt3,lenma3
1877 * Copy of the new m_index from the disk
1878  INTEGER lndex(m_idmax,3)
1879  CHARACTER*80 titld(m_idmax)
1880  LOGICAL glk_exist
1881 *-----------------------------------------------------------
1882  CALL glk_initialize
1883  nouth=m_huni
1884 * Read basic params
1885  READ(nouth,'(6i10)') nvrs3,lenind3,lenma3,lengt3
1886  IF(m_length+lengt3 .GE. m_lenmb) goto 900
1887 * Check version
1888  IF(m_version .NE. nvrs3) WRITE(m_out,*)
1889  $ '+++++WARNING (GLK_hrin): histos produced by older version',nvrs3
1890  IF(m_version .NE. nvrs3) WRITE(6,*)
1891  $ '+++++WARNING (GLK_hrin): histos produced by older version',nvrs3
1892  DO i=1,m_idmax
1893  DO k=1,3
1894  lndex(i,k)=0
1895  ENDDO
1896  ENDDO
1897 * We keep backward compatibility for reading disk files
1898  IF(nvrs3. lt. 130) lenind3 = m_idmax
1899 * Read new index and title from the disk
1900  READ(nouth,'(6i10)') ((lndex(i,k),k=1,3),i=1,lenind3)
1901  READ(nouth,'(a80)') (titld(i), i=1,lenind3)
1902  lenold=m_length
1903 * Append AT ONCE content of new histos at the end of storage m_b
1904  m_length=m_length+lengt3
1905  READ(nouth,'(3d24.16)') (m_b(i),i=lenold+1,m_length)
1906 
1907 * Append ONE BY ONE m_index and m_titlc with new histos
1908  CALL glk_hadres(0,lact)
1909  DO l=1,lenind3
1910  IF(lact .EQ. 0) goto 901
1911  idn= lndex(l,1)
1912  IF(idn .EQ. 0) goto 100
1913 * Identical id's are changed by adding big number = 1000000
1914  10 CONTINUE
1915  IF( glk_exist(idn) ) THEN
1916  idn = idn +1000000*(idn/iabs(idn))
1917  goto 10
1918  ENDIF
1919  m_index(lact,1)=idn
1920  m_index(lact,2)=lndex(l,2)+lenold
1921  m_index(lact,3)=lndex(l,3)
1922  m_titlc(lact) =titld(l)
1923 * Still one small correction in the newly appended histo
1924  istn = m_index(lact,2)
1925  m_b(istn +2) = 9d12 + idn*10 +9d0
1926  lact=lact+1
1927  ENDDO
1928  100 CONTINUE
1929  RETURN
1930 
1931  900 CONTINUE
1932  CALL glk_stop1('++++ GLK_hrin: to litle space, m_LenmB= ',m_lenmb)
1933  901 CONTINUE
1934  CALL glk_stop1('++++ GLK_hrin: to many histos, m_idmax= ',m_idmax)
1935  END
1936 
1937 
1938  SUBROUTINE glk_hrin2(idum1,idum2,idum3)
1939 * **********************************
1940 * New version which has a possibility to
1941 * ADD histograms
1942 * If ID is not existing already then no action is taken
1943 * ***********************************
1944  IMPLICIT NONE
1945  include 'GLK.h'
1946  INTEGER idum1,idum2,idum3
1947 * Copy of the histos from the disk
1948  DOUBLE PRECISION bz(m_lenmb)
1949  INTEGER indez(m_idmax,3)
1950  CHARACTER*80 titlz(m_idmax)
1951  LOGICAL glk_exist
1952  INTEGER nouth,ist3,nchx,ist,ist2,ist3z,nchxz,istz
1953  INTEGER ist2z,lact,lenmaz,lengtz,nvrsz,lenindz,lz,id,i,k
1954 *-------------------------------------------------------------
1955  CALL glk_initialize
1956  nouth=m_huni
1957 * Read basic params
1958  READ(nouth,'(6i10)') nvrsz,lenindz,lenmaz,lengtz
1959 * Check version
1960  IF(m_version .NE. nvrsz) WRITE(m_out,*)
1961  $ '++++WARNING (GLK_hrin2): histos produced by older version',nvrsz
1962  IF(m_version .NE. nvrsz) WRITE(6,*)
1963  $ '++++WARNING (GLK_hrin2): histos produced by older version',nvrsz
1964 
1965 * We keep backward compatibility for reading disk files
1966  IF(nvrsz. lt. 130) lenindz = m_idmax
1967  DO i=1,m_idmax
1968  DO k=1,3
1969  indez(i,k)=0
1970  ENDDO
1971  ENDDO
1972 * Read new m_index, title and bins from the disk
1973  READ(nouth,'(6i10)') ((indez(i,k),k=1,3),i=1,lenindz)
1974  READ(nouth,'(a80)') (titlz(i) , i=1,lenindz)
1975  READ(nouth,'(3d24.16)') (bz(i),i=1,lengtz)
1976 
1977 * Add new histos from disk to existing ones one by one
1978  DO lz=1,lenindz
1979  id= indez(lz,1)
1980  IF(id .EQ. 0) goto 200
1981  IF( .NOT. glk_exist(id)) THEN
1982  CALL glk_retu1('GLK_hrin2: unmached skipped histo ID=', id)
1983  goto 100
1984  ENDIF
1985 * parameters of existing histo
1986  CALL glk_hadres(id,lact)
1987  ist = m_index(lact,2)
1988  ist2 = ist+7
1989  ist3 = ist+11
1990  nchx = m_b(ist2 +1)
1991 * parameters of the histo from the disk
1992  istz = indez(lz,2)
1993  ist2z = istz+7
1994  ist3z = istz+11
1995  nchxz = bz(ist2z +1)
1996  IF(nchx .NE. nchxz) THEN
1997  CALL glk_retu1('GLK_hrin2: non-equal bins, skipped ID=', id)
1998  goto 100
1999  ENDIF
2000 * Add/Merge all additive entries of the two histos
2001 * No of entries in m_index
2002  m_index(lact,3) = m_index(lact,3)+indez(lact,3)
2003 * Overflows, underflows etc.
2004  DO i=1,12
2005  m_b(ist3+i)=m_b(ist3+i) +bz(ist3z+i)
2006  ENDDO
2007 * Except of this one non-additive entry
2008  m_b(ist3+13)=max(m_b(ist3+13),m_b(ist3z+13))
2009 * Regular bin content added now!
2010  DO i= 1, 2*nchx
2011  m_b(ist+m_buf1+i)=m_b(ist+m_buf1+i) +bz(istz+m_buf1+i)
2012  ENDDO
2013  100 CONTINUE
2014  ENDDO
2015  200 CONTINUE
2016  END
2017 
2018  SUBROUTINE glk_hrend(chdum)
2019 * ***********************
2020  IMPLICIT NONE
2021  include 'GLK.h'
2022  CHARACTER*(*) chdum
2023 *---------------------------
2024  CLOSE(m_huni)
2025  END
2026 *======================================================================
2027 * End of histograming procedures
2028 *======================================================================
2029 
2030 
2031 
2032 *======================================================================
2033 * Ploting procedures using LaTeX
2034 *======================================================================
2035 
2036  SUBROUTINE glk_plinitialize(Lint,TeXfile)
2037 *----------------------------------------------------------------------
2038 * Lint =0
2039 * Normal mode, full LaTeX header
2040 * Lint =1
2041 * For TeX file is used in \input, no LaTeX header
2042 * Lint =2
2043 * LaTeX header for one-page plot used as input for postscript
2044 *
2045 * Negative Lint only for debug, big frame around plot is added.
2046 *----------------------------------------------------------------------
2047  IMPLICIT NONE
2048  include 'GLK.h'
2049  SAVE
2050  INTEGER lint,noufig
2051  CHARACTER*60 texfile
2052 *---------------------------------
2053 * Initialize GLK_Plot
2054  CALL glk_plint(lint) ! Define header style
2055  noufig=11 !
2056  OPEN(noufig,file=texfile) ! Open LaTeX file
2057  CALL glk_plcap(noufig) ! Initialize GLK_Plot
2058  END
2059 
2060  SUBROUTINE glk_plend
2061 * ********************
2062  IMPLICIT NONE
2063  include 'GLK.h'
2064  SAVE
2065 *---------------------------------------------------
2066 * Note that TeX file is used in \input then you may not want
2067 * to have header and \end{document}
2068  IF( abs(m_lint) .NE. 1) THEN
2069  WRITE(m_ltx,'(2A)') m_bs,'end{document}'
2070  ENDIF
2071  CLOSE(m_ltx)
2072  END
2073 
2074  SUBROUTINE glk_plint(Lint)
2075 * **************************
2076  IMPLICIT NONE
2077  include 'GLK.h'
2078  SAVE
2079  INTEGER lint
2080 *---------------------------------
2081  m_lint = lint
2082  END
2083 
2084  SUBROUTINE glk_plcap(LtxUnit)
2085 * ****************************
2086  IMPLICIT NONE
2087  include 'GLK.h'
2088  INTEGER ltxunit,i,k
2089 *---------
2090  CALL glk_initialize
2091  m_keytit= 0
2092  DO i=1,m_titlen
2093  DO k=1,80
2094  m_titch(i)(k:k)=' '
2095  ENDDO
2096  ENDDO
2097 *---------
2098  m_tline = 1
2099  m_ltx=iabs(ltxunit)
2100 
2101  IF( abs(m_lint) .EQ. 0) THEN
2102 * Normal mode, no colors!!!
2103  WRITE(m_ltx,'(A,A)') m_bs,'documentclass[12pt]{article}'
2104 *!! WRITE(m_ltx,'(A,A)') m_BS,'usepackage{html}'
2105  WRITE(m_ltx,'(A,A)') m_bs,'textwidth = 16cm'
2106  WRITE(m_ltx,'(A,A)') m_bs,'textheight = 18cm'
2107  WRITE(m_ltx,'(A,A)') m_bs,'begin{document}'
2108  WRITE(m_ltx,'(A)') ' '
2109  ELSEIF( abs(m_lint) .EQ. 1) THEN
2110 * For TeX file is used in \input
2111  WRITE(m_ltx,'(A)') ' '
2112  ELSEIF( abs(m_lint) .EQ. 2) THEN
2113 * For one-page plot being input for postrscript
2114 *!! WRITE(m_ltx,'(A,A)') m_BS,'documentclass[12pt,html]{article}'
2115 *!! WRITE(m_ltx,'(A,A)') m_BS,'documentclass[12pt,dvips]{seminar}' !<-for colors!!!
2116  WRITE(m_ltx,'(A,A)') m_bs,'documentclass[12pt,dvips]{article}'
2117  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{amsmath}'
2118  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{amssymb}'
2119 *!! WRITE(m_ltx,'(A,A)') m_BS,'usepackage{html}'
2120  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{epsfig}'
2121  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{epic}'
2122  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{eepic}'
2123  WRITE(m_ltx,'(A,A)') m_bs,'usepackage{color}' !<-for colors!!!
2124 *!! WRITE(m_ltx,'(A,A)') m_BS,'hoffset = -1in'
2125 *!! WRITE(m_ltx,'(A,A)') m_BS,'voffset = -1in'
2126 *!! WRITE(m_ltx,'(A,A)') m_BS,'textwidth = 16cm'
2127 *!! WRITE(m_ltx,'(A,A)') m_BS,'textheight = 16cm'
2128 *!! WRITE(m_ltx,'(A,A)') m_BS,'oddsidemargin = 0cm'
2129 *!! WRITE(m_ltx,'(A,A)') m_BS,'topmargin = 0cm'
2130 *!! WRITE(m_ltx,'(A,A)') m_BS,'headheight = 0cm'
2131 *!! WRITE(m_ltx,'(A,A)') m_BS,'headsep = 0cm'
2132  WRITE(m_ltx,'(A,A)') m_bs,'begin{document}'
2133  WRITE(m_ltx,'(A,A)') m_bs,'pagestyle{empty}'
2134  WRITE(m_ltx,'(A)') ' '
2135  ELSE
2136  CALL glk_stop1('+++STOP in GLK_PlInt, wrong m_lint=',m_lint)
2137  ENDIF
2138  END
2139 
2140 
2141  SUBROUTINE glk_plot(id,ch1,ch2,kdum)
2142 * ************************************
2143  IMPLICIT NONE
2144  include 'GLK.h'
2145  CHARACTER ch1,ch2,chr
2146  CHARACTER*80 title
2147  INTEGER id,kdum
2148  DOUBLE PRECISION yy(m_maxnb),yer(m_maxnb)
2149  LOGICAL glk_exist
2150  INTEGER idum,kax,kay,ioplog,iopsla,ioperb,iopsc1,iopsc2
2151  INTEGER ker,nchx
2152  DOUBLE PRECISION xl,xu,dxl,dxu,yl,yu
2153 *--------------------------------------------
2154  DATA chr /' '/
2155 * RETURN if histo non-existing
2156  IF(.NOT.glk_exist(id)) goto 900
2157 * ...unpack histogram
2158  CALL glk_unpak(id,yy ,' ',idum)
2159  CALL glk_unpak(id,yer,'ERRO',idum)
2160  CALL glk_hinbo1(id,title,nchx,dxl,dxu)
2161  xl = dxl
2162  xu = dxu
2163  CALL glk_range1(id,yl,yu)
2164  kax=1200
2165  kay=1200
2166  IF(ch1 .EQ. 'S') THEN
2167 * ...superimpose plot
2168  backspace(m_ltx)
2169  backspace(m_ltx)
2170  ELSE
2171 * ...new frame only
2172  chr=ch1
2173  CALL glk_plfram1(id,kax,kay)
2174  ENDIF
2175  WRITE(m_ltx,'(A)') '%========== next plot (line) =========='
2176  WRITE(m_ltx,'(A,I10)') '%==== HISTOGRAM ID=',id
2177  WRITE(m_ltx,'(A,A70 )')'% ',title
2178 *...cont. line for functions
2179  CALL glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
2180  ker = ioperb-1
2181  IF (iopsla .EQ. 2) chr='C'
2182 *...suppress GLK_PLOT assignments
2183  IF (ch2 .EQ. 'B') chr=' '
2184  IF (ch2 .EQ. '*') chr='*'
2185  IF (ch2 .EQ. 'C') chr='C'
2186 *...various types of lines
2187  IF (chr .EQ. ' ') THEN
2188 *...contour line used for histogram
2189  CALL glk_plhist(kax,kay,nchx,yl,yu,yy,ker,yer)
2190  ELSE IF(chr .EQ. '*') THEN
2191 *...marks in the midle of the bin
2192  CALL glk_plhis2(kax,kay,nchx,yl,yu,yy,ker,yer)
2193  ELSE IF(chr .EQ. 'C') THEN
2194 *...slanted (dotted) line in plotting non-MC functions
2195  CALL glk_plcirc(kax,kay,nchx,yl,yu,yy)
2196  ENDIF
2197 *------------------------------!
2198 * Ending
2199 *------------------------------!
2200  WRITE(m_ltx,'(2A)') m_bs,'end{picture} % close entire picture '
2201  WRITE(m_ltx,'(2A)') m_bs,'end{figure}'
2202 
2203  RETURN
2204  900 CALL glk_retu1('+++GLK_PLOT: Nonexistig histo, skipped, id=' ,id)
2205  END
2206 
2207  SUBROUTINE glk_plfram1(ID,kax,kay)
2208 * **********************************
2209  IMPLICIT NONE
2210  include 'GLK.h'
2211  INTEGER id,kax,kay
2212  CHARACTER*80 title
2213  DOUBLE PRECISION tipsy(20),tipsx(20)
2214  DOUBLE PRECISION xl,dxl,xu,dxu
2215  INTEGER ntipy,ntipx,nchx,icont
2216  DOUBLE PRECISION yu,yl
2217  DATA icont/0/
2218 *----------------
2219  icont=icont+1
2220  CALL glk_hinbo1(id,title,nchx,dxl,dxu)
2221  xl = dxl
2222  xu = dxu
2223  CALL glk_range1(id,yl,yu)
2224 
2225  IF(icont .GT. 1) WRITE(m_ltx,'(2A)') m_bs,'newpage'
2226 *------------------------------!
2227 * Header
2228 *------------------------------!
2229  WRITE(m_ltx,'(A)') ' '
2230  WRITE(m_ltx,'(A)') ' '
2231  WRITE(m_ltx,'(A)') '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2232  $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
2233  WRITE(m_ltx,'(A)') '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2234  $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
2235  WRITE(m_ltx,'(2A)') m_bs,'begin{figure}[!ht]'
2236  WRITE(m_ltx,'(2A)') m_bs,'centering'
2237 *------------------------------!
2238 * General Caption
2239 *------------------------------!
2240  WRITE(m_ltx,'(4A)') m_bs,'caption{',m_bs,'small'
2241  IF(m_keytit.EQ.0) THEN
2242  WRITE(m_ltx,'(A)') title
2243  ELSE
2244  WRITE(m_ltx,'(A)') m_titch(1)
2245  ENDIF
2246  WRITE(m_ltx,'(A)') '}'
2247 *------------------------------!
2248 * Frames and labels
2249 *------------------------------!
2250  WRITE(m_ltx,'(A)') '% =========== big frame, title etc. ======='
2251  WRITE(m_ltx,'(4A)') m_bs,'setlength{',m_bs,'unitlength}{0.1mm}'
2252  WRITE(m_ltx,'(2A)') m_bs,'begin{picture}(1600,1500)'
2253  WRITE(m_ltx,'(4A)')
2254  $ m_bs,'put(0,0){',m_bs,'framebox(1600,1500){ }}'
2255  WRITE(m_ltx,'(A)') '% =========== small frame, labeled axis ==='
2256  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2257  $ m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
2258  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2259  $ m_bs,'put(0,0){',m_bs,'framebox( ',kax,',',kay,'){ }}'
2260  WRITE(m_ltx,'(A)') '% =========== x and y axis ================'
2261  CALL glk_saxisx(kax,xl,xu,ntipx,tipsx)
2262  CALL glk_saxisy(kay,yl,yu,ntipy,tipsy)
2263  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}'
2264  $ ,'% end of plotting labeled axis'
2265  END
2266 
2267  SUBROUTINE glk_saxisx(kay,YL,YU,NLT,TIPSY)
2268 * ***************************************
2269 * plotting x-axis with long and short tips
2270  IMPLICIT NONE
2271  include 'GLK.h'
2272  INTEGER kay,nlt
2273  DOUBLE PRECISION yl,yu,tipsy(20)
2274 *
2275  INTEGER ly,jy,n,nts,k,lex
2276  DOUBLE PRECISION dy,pds,scmx,p0s,ddys,yy0l,ddyl,pdl,p0l,yy0s
2277 *---------------------------------------------------
2278  dy= abs(yu-yl)
2279  ly = nint( log10(dy) -0.4999999d0 )
2280  jy = nint(dy/10d0**ly)
2281  ddyl = dy*10d0**(-ly)
2282  IF( jy .EQ. 1) ddyl = 10d0**ly*0.25d0
2283  IF( jy .GE. 2.AND.jy .LE. 3) ddyl = 10d0**ly*0.5d0
2284  IF( jy .GE. 4.AND.jy .LE. 6) ddyl = 10d0**ly*1.0d0
2285  IF( jy .GE. 7) ddyl = 10d0**ly*2.0d0
2286  WRITE(m_ltx,'(A)') '% .......GLK_SAxisX........ '
2287  WRITE(m_ltx,'(A,I4)') '% JY= ',jy
2288 *-------
2289  nlt = int(dy/ddyl)
2290  nlt = max0(min0(nlt,20),1)+1
2291  yy0l = nint(yl/ddyl+0.5d0)*ddyl
2292  ddys = ddyl/10d0
2293  yy0s = nint(yl/ddys+0.4999999d0)*ddys
2294  p0l = kay*(yy0l-yl)/(yu-yl)
2295  pdl = kay*ddyl/(yu-yl)
2296  p0s = kay*(yy0s-yl)/(yu-yl)
2297  pds = kay*ddys/(yu-yl)
2298  nlt = int(abs(yu-yy0l)/ddyl+0.0000001d0)+1
2299  nts = int(abs(yu-yy0s)/ddys+0.0000001d0)+1
2300  DO 41 n=1,nlt
2301  tipsy(n) =yy0l+ ddyl*(n-1)
2302  41 CONTINUE
2303  WRITE(m_ltx,1000)
2304  $ m_bs,'multiput(' ,p0l, ',0)(' ,pdl, ',0){' ,nlt, '}{',
2305  $ m_bs,'line(0,1){25}}',
2306  $ m_bs,'multiput(' ,p0s, ',0)(' ,pds, ',0){' ,nts, '}{',
2307  $ m_bs,'line(0,1){10}}'
2308  WRITE(m_ltx,1001)
2309  $ m_bs,'multiput(' ,p0l, ',' ,kay, ')(' ,pdl, ',0){' ,nlt,
2310  $ '}{' ,m_bs, 'line(0,-1){25}}',
2311  $ m_bs,'multiput(' ,p0s, ',' ,kay, ')(' ,pds, ',0){' ,nts,
2312  $ '}{' ,m_bs, 'line(0,-1){10}}'
2313  1000 FORMAT(2a,f8.2,a,f8.2,a,i4,3a)
2314  1001 FORMAT(2a,f8.2,a,i4,a,f8.2,a,i4,3a)
2315 * ...labeling of axis
2316  scmx = dmax1(dabs(yl),dabs(yu))
2317  lex = nint( log10(scmx) -0.50001)
2318  DO 45 n=1,nlt
2319  k = nint(kay*(tipsy(n)-yl)/(yu-yl))
2320  IF(lex .LT. 2.AND.lex .GT. -1) THEN
2321 * ...without exponent
2322  WRITE(m_ltx,'(2A,I4,5A,F8.3,A)')
2323  $ m_bs,'put(',k,',-25){',m_bs,'makebox(0,0)[t]{',m_bs,'large $ ',
2324  $ tipsy(n), ' $}}'
2325  ELSE
2326 * ...with exponent
2327  WRITE(m_ltx,'(2A,I4,5A,F8.3,2A,I4,A)')
2328  $ m_bs,'put(' ,k, ',-25){',m_bs,'makebox(0,0)[t]{',m_bs,'large $ ',
2329  $ tipsy(n)/(10d0**lex),m_bs,'cdot 10^{',lex,'} $}}'
2330  ENDIF
2331  45 CONTINUE
2332  END
2333 
2334  SUBROUTINE glk_saxisy(kay,yl,yu,nlt,tipsy)
2335 * ******************************************
2336 * plotting y-axis with long and short tips
2337  IMPLICIT NONE
2338  include 'GLK.h'
2339  INTEGER kay,nlt
2340  DOUBLE PRECISION yl,yu,tipsy(20)
2341 *
2342  DOUBLE PRECISION dy,ddyl,z0l,scmx,yy0s,ddys,yy0l,p0l,pds,p0s,pdl
2343  INTEGER ly,jy,n,nts,k,lex
2344 *---------------------------------------------------
2345  dy= abs(yu-yl)
2346  ly = nint( log10(dy) -0.49999999d0 )
2347  jy = nint(dy/10d0**ly)
2348  ddyl = dy*10d0**(-ly)
2349  IF( jy .EQ. 1) ddyl = 10d0**ly*0.25d0
2350  IF( jy .GE. 2.AND.jy .LE. 3) ddyl = 10d0**ly*0.5d0
2351  IF( jy .GE. 4.AND.jy .LE. 6) ddyl = 10d0**ly*1.0d0
2352  IF( jy .GE. 7) ddyl = 10d0**ly*2.0d0
2353  WRITE(m_ltx,'(A)') '% .......GLK_SAxisY........ '
2354  WRITE(m_ltx,'(A,I4)') '% JY= ',jy
2355 *-------
2356  nlt = int(dy/ddyl)
2357  nlt = max0(min0(nlt,20),1)+1
2358  yy0l = nint(yl/ddyl+0.4999999d0)*ddyl
2359  ddys = ddyl/10d0
2360  yy0s = nint(yl/ddys+0.5d0)*ddys
2361  p0l = kay*(yy0l-yl)/(yu-yl)
2362  pdl = kay*ddyl/(yu-yl)
2363  p0s = kay*(yy0s-yl)/(yu-yl)
2364  pds = kay*ddys/(yu-yl)
2365  nlt= int(abs(yu-yy0l)/ddyl+0.0000001d0) +1
2366  nts= int(abs(yu-yy0s)/ddys+0.0000001d0) +1
2367  DO 41 n=1,nlt
2368  tipsy(n) =yy0l+ ddyl*(n-1)
2369  41 CONTINUE
2370 * plotting tics on vertical axis
2371  WRITE(m_ltx,1000)
2372  $ m_bs,'multiput(0,' ,p0l, ')(0,' ,pdl ,'){' ,nlt, '}{',
2373  $ m_bs,'line(1,0){25}}',
2374  $ m_bs,'multiput(0,' ,p0s, ')(0,' ,pds, '){' ,nts, '}{',
2375  $ m_bs,'line(1,0){10}}'
2376  WRITE(m_ltx,1001)
2377  $ m_bs,'multiput(' ,kay, ',' ,p0l, ')(0,' ,pdl, '){' ,nlt,
2378  $ '}{',m_bs,'line(-1,0){25}}',
2379  $ m_bs,'multiput(' ,kay, ',' ,p0s, ')(0,' ,pds, '){' ,nts,
2380  $ '}{',m_bs,'line(-1,0){10}}'
2381  1000 FORMAT(2a,f8.2,a,f8.2,a,i4,3a)
2382  1001 FORMAT(2a,i4,a,f8.2,a,f8.2,a,i4,3a)
2383 * ...Zero line if necessary
2384  z0l = kay*(-yl)/(yu-yl)
2385  IF(z0l .GT. 0d0.AND.z0l .LT. float(kay))
2386  $ WRITE(m_ltx,'(2A,F8.2,3A,I4,A)')
2387  $ m_bs,'put(0,' ,z0l, '){',m_bs,'line(1,0){' ,kay, '}}'
2388 * ...labeling of axis
2389  scmx = dmax1(dabs(yl),dabs(yu))
2390  lex = nint( log10(scmx) -0.50001d0)
2391  DO 45 n=1,nlt
2392  k = nint(kay*(tipsy(n)-yl)/(yu-yl))
2393  IF(lex .LT. 2.AND.lex .GT. -1) THEN
2394 * ...without exponent
2395  WRITE(m_ltx,'(2A,I4,5A,F8.3,A)')
2396  $ m_bs,'put(-25,' ,k, '){',m_bs,'makebox(0,0)[r]{',
2397  $ m_bs,'large $ ' ,tipsy(n), ' $}}'
2398  ELSE
2399 * ...with exponent
2400  WRITE(m_ltx,'(2A,I4,5A,F8.3,2A,I4,A)')
2401  $ m_bs,'put(-25,' ,k, '){',m_bs,'makebox(0,0)[r]{',
2402  $ m_bs,'large $ '
2403  $ ,tipsy(n)/(10d0**lex), m_bs,'cdot 10^{' ,lex, '} $}}'
2404  ENDIF
2405  45 CONTINUE
2406  END
2407 
2408  SUBROUTINE glk_plhist(kax,kay,nchx,yl,yu,yy,ker,yer)
2409 * ************************************************
2410 * plotting contour line for histogram
2411 * ***********************
2412  IMPLICIT NONE
2413  include 'GLK.h'
2414  INTEGER kax,kay,nchx,ker
2415  DOUBLE PRECISION yl,yu,yy(*),yer(*)
2416  CHARACTER*80 fmt1
2417 *
2418  INTEGER ix0,ix2,idx,ie,ierr,idy,ib,iy0,iy1,ix1
2419 *---------------------------------------------------
2420  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2421  $ m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
2422  WRITE(m_ltx,'(A)') '% ========== plotting primitives =========='
2423 *...various types of line
2424  IF(m_tline .EQ. 1) THEN
2425  WRITE(m_ltx,'(2A)') m_bs,'thicklines '
2426  ELSE
2427  WRITE(m_ltx,'(2A)') m_bs,'thinlines '
2428  ENDIF
2429 *...short macros for vertical/horizontal straight lines
2430  WRITE(m_ltx,'(8A)')
2431  $ m_bs,'newcommand{',m_bs,'x}[3]{',m_bs,'put(#1,#2){',
2432  $ m_bs,'line(1,0){#3}}}'
2433  WRITE(m_ltx,'(8A)')
2434  $ m_bs,'newcommand{',m_bs,'y}[3]{',m_bs,'put(#1,#2){',
2435  $ m_bs,'line(0,1){#3}}}'
2436  WRITE(m_ltx,'(8A)')
2437  $ m_bs,'newcommand{',m_bs,'z}[3]{',m_bs,'put(#1,#2){',
2438  $ m_bs,'line(0,-1){#3}}}'
2439 * error bars
2440  WRITE(m_ltx,'(8A)')
2441  $ m_bs,'newcommand{',m_bs,'e}[3]{',
2442  $ m_bs,'put(#1,#2){',m_bs,'line(0,1){#3}}}'
2443  ix0=0
2444  iy0=0
2445  DO 100 ib=1,nchx
2446  ix1 = nint(kax*(ib-0.00001)/nchx) !ib=7
2447  iy1 = nint(kay*(yy(ib)-yl)/(yu-yl)) !iy1=775,while ix0=168,iy0=770
2448  idy = iy1-iy0
2449  idx = ix1-ix0
2450  fmt1 = '(2(2A,I4,A,I4,A,I4,A))'
2451  IF( idy .GE. 0) THEN
2452  IF(iy1 .GE. 0.AND.iy1 .LE. kay)
2453  $ WRITE(m_ltx,fmt1) m_bs,'y{',ix0,'}{',iy0,'}{',idy,'}',
2454  $ m_bs,'x{',ix0,'}{',iy1,'}{',idx,'}'
2455  ELSE
2456  IF(iy1 .GE. 0.AND.iy1 .LE. kay)
2457  $ WRITE(m_ltx,fmt1) m_bs,'z{',ix0,'}{',iy0,'}{',-idy,'}',
2458  $ m_bs,'x{',ix0,'}{',iy1,'}{',idx,'}'
2459  ENDIF
2460  ix0=ix1
2461  iy0=iy1
2462  IF(ker .EQ. 1) THEN
2463  ix2 = nint(kax*(ib-0.5000d0)/nchx)
2464  ierr = nint(kay*((yy(ib)-yer(ib))-yl)/(yu-yl))
2465  ie = nint(kay*yer(ib)/(yu-yl))
2466  IF(iy1 .GE. 0.AND.iy1 .LE. kay.and.abs(ierr) .LE. 9999
2467  $ .and.2*ie .LE. 9999) WRITE(m_ltx,8000) m_bs,ix2,ierr,ie*2
2468  ENDIF
2469  100 CONTINUE
2470 8000 FORMAT(4(a1,2he{,i4,2h}{,i5,2h}{,i4,1h}:1x ))
2471  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}',
2472  $ ' % end of plotting histogram'
2473 * change line-style
2474  m_tline= m_tline+1
2475  IF(m_tline .GT. 2) m_tline=1
2476  END
2477 
2478  SUBROUTINE glk_plhis2(kax,kay,nchx,yl,yu,yy,ker,yer)
2479 * ************************************************
2480 * marks in the midle of the bin
2481 * **********************************
2482  IMPLICIT NONE
2483  include 'GLK.h'
2484  DOUBLE PRECISION yl,yu,yy(*),yer(*)
2485  INTEGER kax,kay,nchx,ker
2486 *
2487  INTEGER iy1,ierr,ie,ix1,irad1,irad2,ib
2488 *---------------------------------------------------
2489 
2490  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2491  $ m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
2492  WRITE(m_ltx,'(A)') '% ========== plotting primitives =========='
2493 *...various types of mark
2494  irad1= 6
2495  irad2=10
2496  IF(m_tline .EQ. 1) THEN
2497 * small filled circle
2498  WRITE(m_ltx,'(8A,I3,A)')
2499  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2500  $ m_bs,'put(#1,#2){',m_bs,'circle*{',irad1,'}}}'
2501  ELSEIF(m_tline .EQ. 2) THEN
2502 * small open circle
2503  WRITE(m_ltx,'(8A,I3,A)')
2504  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2505  $ m_bs,'put(#1,#2){',m_bs,'circle{',irad1,'}}}'
2506  ELSEIF(m_tline .EQ. 3) THEN
2507 * big filled circle
2508  WRITE(m_ltx,'(8A,I3,A)')
2509  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2510  $ m_bs,'put(#1,#2){',m_bs,'circle*{',irad2,'}}}'
2511  ELSEIF(m_tline .EQ. 4) THEN
2512 * big open circle
2513  WRITE(m_ltx,'(8A,I3,A)')
2514  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2515  $ m_bs,'put(#1,#2){',m_bs,'circle{',irad2,'}}}'
2516 * Other symbols
2517  ELSEIF(m_tline .EQ. 5) THEN
2518  WRITE(m_ltx,'(10A)')
2519  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2520  $ m_bs,'put(#1,#2){',m_bs,'makebox(0,0){$',m_bs,'diamond$}}}'
2521  ELSE
2522  WRITE(m_ltx,'(10A)')
2523  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2524  $ m_bs,'put(#1,#2){',m_bs,'makebox(0,0){$',m_bs,'star$}}}'
2525  ENDIF
2526 * error bars
2527  WRITE(m_ltx,'(8A)')
2528  $ m_bs,'newcommand{',m_bs,'E}[3]{',
2529  $ m_bs,'put(#1,#2){',m_bs,'line(0,1){#3}}}'
2530  DO 100 ib=1,nchx
2531  ix1 = nint(kax*(ib-0.5000d0)/nchx)
2532  iy1 = nint(kay*(yy(ib)-yl)/(yu-yl))
2533  IF(iy1 .GE. 0.AND.iy1 .LE. kay) WRITE(m_ltx,7000) m_bs,ix1,iy1
2534  IF(ker .EQ. 1) THEN
2535  ierr = nint(kay*((yy(ib)-yer(ib))-yl)/(yu-yl))
2536  ie = nint(kay*yer(ib)/(yu-yl))
2537  IF(iy1 .GE. 0.AND.iy1 .LE. kay.and.abs(ierr) .LE. 9999
2538  $ .and.2*ie .LE. 9999) WRITE(m_ltx,8000) m_bs,ix1,ierr,ie*2
2539  ENDIF
2540  100 CONTINUE
2541 7000 FORMAT(4(a1,2hr{,i4,2h}{,i4,1h}:1x ))
2542 8000 FORMAT(4(a1,2he{,i4,2h}{,i5,2h}{,i4,1h}:1x ))
2543  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}',
2544  $ ' % end of plotting histogram'
2545 * change line-style
2546  m_tline= m_tline+1
2547  IF(m_tline .GT. 6) m_tline=1
2548  END
2549 
2550  SUBROUTINE glk_plcirc(kax,kay,nchx,yl,yu,yy)
2551 * ****************************************
2552 * plots equidistant points, four-point interpolation,
2553  IMPLICIT NONE
2554  include 'GLK.h'
2555  INTEGER kax,kay,nchx
2556  DOUBLE PRECISION yl,yu,yy(*)
2557 *
2558  INTEGER ix(m_maxnb),iy(m_maxnb)
2559  DOUBLE PRECISION ai0,dx,aj0,ds,facy,aj,ai
2560  INTEGER ipnt,i,iter,ipoin,irad1,irad2
2561  DOUBLE PRECISION glk_aprof
2562 *---------------------------------------------------
2563 
2564 * ...various types of line
2565 * ...distance between points is DS, radius of a point is IRAD
2566  irad2=6
2567  irad1=3
2568 * .............
2569  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2570  $ m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
2571  WRITE(m_ltx,'(A)') '% ========== plotting primitives =========='
2572  IF(m_tline .EQ. 1) THEN
2573 * small filled circle
2574  ds = 10
2575  WRITE(m_ltx,'(8A,I3,A)')
2576  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2577  $ m_bs,'put(#1,#2){',m_bs,'circle*{',irad1,'}}}'
2578  ELSEIF(m_tline .EQ. 2) THEN
2579 * small open circle
2580  ds = 10
2581  WRITE(m_ltx,'(8A,I3,A)')
2582  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2583  $ m_bs,'put(#1,#2){',m_bs,'circle{',irad1,'}}}'
2584  ELSEIF(m_tline .EQ. 3) THEN
2585 * big filled circle
2586  ds = 20
2587  WRITE(m_ltx,'(8A,I3,A)')
2588  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2589  $ m_bs,'put(#1,#2){',m_bs,'circle*{',irad2,'}}}'
2590  ELSEIF(m_tline .EQ. 4) THEN
2591 * big open circle
2592  ds = 20
2593  WRITE(m_ltx,'(8A,I3,A)')
2594  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2595  $ m_bs,'put(#1,#2){',m_bs,'circle{',irad2,'}}}'
2596 * Other symbols
2597  ELSEIF(m_tline .EQ. 5) THEN
2598  ds = 20
2599  WRITE(m_ltx,'(10A)')
2600  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2601  $ m_bs,'put(#1,#2){',m_bs,'makebox(0,0){$',m_bs,'diamond$}}}'
2602  ELSE
2603  ds = 20
2604  WRITE(m_ltx,'(10A)')
2605  $ m_bs,'newcommand{',m_bs,'R}[2]{',
2606  $ m_bs,'put(#1,#2){',m_bs,'makebox(0,0){$',m_bs,'star$}}}'
2607  ENDIF
2608  facy = kay/(yu-yl)
2609 * plot first point
2610  ai = 0.
2611  aj = (glk_aprof( (ai/kax)*nchx+0.5d0, nchx, yy) -yl)*facy
2612  ipnt =1
2613  ix(ipnt) = int(ai)
2614  iy(ipnt) = int(aj)
2615  dx = ds
2616  ai0 = ai
2617  aj0 = aj
2618 * plot next points
2619  DO 100 ipoin=2,3000
2620 * iteration to get (approximately) equal distance among ploted points
2621  DO 50 iter=1,3
2622  ai = ai0+dx
2623  aj = (glk_aprof( (ai/kax)*nchx+0.5d0, nchx, yy) -yl)*facy
2624  dx = dx *ds/sqrt(dx**2 + (aj-aj0)**2)
2625  50 CONTINUE
2626  IF(int(aj) .GE. 0.AND.int(aj) .LE. kay.AND.int(ai) .LE. kax) THEN
2627  ipnt = ipnt+1
2628  ix(ipnt) = int(ai)
2629  iy(ipnt) = int(aj)
2630  ENDIF
2631  ai0 = ai
2632  aj0 = aj
2633  IF(int(ai) .GT. kax) goto 101
2634  100 CONTINUE
2635  101 CONTINUE
2636  WRITE(m_ltx,7000) (m_bs,ix(i),iy(i), i=1,ipnt)
2637 7000 FORMAT(4(a1,2hr{,i4,2h}{,i4,1h}:1x ))
2638  WRITE(m_ltx,'(2A)') m_bs,'end{picture}} % end of plotting line'
2639 * change line-style
2640  m_tline= m_tline+1
2641  IF(m_tline .GT. 2) m_tline=1
2642  END
2643 
2644  DOUBLE PRECISION FUNCTION glk_aprof(px,nch,yy)
2645 * ************************************************
2646 * PX is a continuous extension of the m_index in array YY
2647  IMPLICIT NONE
2648  INTEGER nch,ip
2649  DOUBLE PRECISION px,yy(*),x,p
2650 *-----------------------------------------------------
2651  x=px
2652  IF(x .LT. 0.0.OR.x .GT. float(nch+1)) THEN
2653  glk_aprof= -1e-20
2654  RETURN
2655  ENDIF
2656  ip=int(x)
2657  IF(ip .LT. 2) ip=2
2658  IF(ip .GT. nch-2) ip=nch-2
2659  p=x-ip
2660  glk_aprof =
2661  $ -(1./6.)*p*(p-1)*(p-2) *yy(ip-1)
2662  $ +(1./2.)*(p*p-1)*(p-2) *yy(ip )
2663  $ -(1./2.)*p*(p+1)*(p-2) *yy(ip+1)
2664  $ +(1./6.)*p*(p*p-1) *yy(ip+2)
2665  END
2666 
2667  SUBROUTINE glk_pltitle(title)
2668 * *****************************
2669  IMPLICIT NONE
2670  include 'GLK.h'
2671  SAVE
2672  CHARACTER*80 title
2673 *----------------------------------------
2674  m_keytit=1
2675  CALL glk_copch(title,m_titch(1))
2676  END
2677 
2678  SUBROUTINE glk_plcapt(lines)
2679 * ****************************
2680 * This routine defines caption and should be called
2681 * before CALL GLK_Plot2, GLK_PlTable or bpltab2
2682 * The matrix CHARACTER*80 lines containes text of the caption ended
2683 * with the last line '% end-of-caption'
2684  IMPLICIT NONE
2685  CHARACTER*80 lines(*)
2686  include 'GLK.h'
2687  SAVE
2688  INTEGER i
2689 *----------------------------------
2690  m_keytit=0
2691  DO i=1,m_titlen
2692  m_titch(i)=lines(i)
2693  m_keytit= m_keytit+1
2694  IF(lines(i) .EQ. '% end-of-caption' ) goto 100
2695  ENDDO
2696  CALL glk_retu1(' WARNING from GLK_PlCapt: to many lines =',m_titlen)
2697  100 CONTINUE
2698  END
2699 
2700  SUBROUTINE glk_pllabel(lines)
2701 * *****************************
2702 * This should be envoked after 'CALL GLK_Plot2'
2703 * to add lines of TeX to a given plot
2704 * ***********************************
2705  IMPLICIT NONE
2706  CHARACTER*80 lines(*)
2707  include 'GLK.h'
2708  SAVE
2709  INTEGER i
2710 *----------------------------------
2711  m_keytit=0
2712  DO i=1,m_titlen
2713  m_titch(i)=lines(i)
2714  m_keytit= m_keytit+1
2715  IF(lines(i) .EQ. '% end-of-label' ) goto 100
2716  ENDDO
2717  CALL glk_retu1(' WARNING from GLK_PlLabel: to many lines =',m_titlen)
2718  100 CONTINUE
2719 *------------------------------!
2720 * erase Ending !
2721 *------------------------------!
2722  backspace(m_ltx)
2723  backspace(m_ltx)
2724 *
2725  DO i=1,m_keytit
2726  WRITE(m_ltx,'(A)') m_titch(i)
2727  ENDDO
2728 *------------------------------!
2729 * restore Ending !
2730 *------------------------------!
2731  WRITE(m_ltx,'(2A)') m_bs,'end{picture} % close entire picture '
2732  IF(abs(m_lint) .EQ. 2) THEN
2733  WRITE(m_ltx,'(A)') '%====== end of GLK_PlLabel =========='
2734  ELSE
2735  WRITE(m_ltx,'(2A)') m_bs,'end{figure}'
2736  ENDIF
2737  END
2738 
2739 
2740  SUBROUTINE glk_plot2(id,ch1,ch2,chmark,chxfmt,chyfmt)
2741 * *****************************************************
2742 * The new, more user-friendly, version of older GLK_Plot
2743 * INPUT:
2744 * ID histogram identifier
2745 * ch1 = ' ' normal new plot
2746 * = 'S' impose new plot on previous one
2747 * ch2 = ' ' ploting line default, contour
2748 * = '*' error bars in midle of the bin
2749 * = 'R' error bars at Right edge of the bin
2750 * = 'L' error bars at Left edge of the bin
2751 * = 'C' slanted continuous smooth line
2752 * chmark = TeX symbol for ploting points
2753 * chxfmt = format (string) for labeling x-axis
2754 * chyfmt = format (string) for labeling y-axis
2755 * Furthermore:
2756 * Captions are defined by means of
2757 * CALL GLK_PlCapt(capt) before CALL GLK_Plot2
2758 * where CHARACTER*80 capt(50) is content of
2759 * caption, line by line, see also comments in GLK_PlCapt routine.
2760 * Additional text as a TeX source text can be appended by means of
2761 * CALL GLK_PlLabel(lines) after CALL GLK_Plot2
2762 * where CHARACTER*80 lines(50) is the TeX add-on.
2763 * This is to be used to decorate plot with
2764 * any kind marks, special labels and text on the plot.
2765 *
2766 * ************************************
2767  IMPLICIT NONE
2768  INTEGER id
2769  CHARACTER ch1,ch2,chmark*(*)
2770  CHARACTER*8 chxfmt,chyfmt
2771  include 'GLK.h'
2772  SAVE
2773  DOUBLE PRECISION yy(m_maxnb),yer(m_maxnb)
2774  CHARACTER*80 title
2775 *---------------------------------------------------------------------
2776  LOGICAL glk_exist
2777  INTEGER kax,kay,incr,ker,nchx
2778  INTEGER iopsla,ioplog,ioperb,iopsc1,iopsc2,idum
2779  DOUBLE PRECISION dxl,dxu,xu,xl,yu,yl
2780  CHARACTER chr
2781  DATA chr /' '/
2782 * TeX Names of the error-bar command and of the point-mark command
2783  CHARACTER*1 chre, chrp1
2784  parameter( chre = 'E', chrp1= 'R' )
2785  CHARACTER*2 chrp
2786 * TeX Name of the point-mark command
2787  CHARACTER*1 chrx(12)
2788  DATA chrx /'a','b','c','d','f','g','h','i','j','k','l','m'/
2789 *---------------------------------------------------------------------
2790 * RETURN if histo non-existing
2791  IF(.NOT.glk_exist(id)) goto 900
2792 * ...unpack histogram
2793  CALL glk_unpak(id,yy ,' ',idum)
2794  CALL glk_unpak(id,yer,'ERRO',idum)
2795  CALL glk_hinbo1(id,title,nchx,dxl,dxu)
2796 * Header
2797  kax=1200
2798  kay=1200
2799  IF(ch1 .EQ. 'S') THEN
2800 * Superimpose plot
2801  incr=incr+1
2802  backspace(m_ltx)
2803  backspace(m_ltx)
2804  ELSE
2805 * New frame only
2806  incr=1
2807  chr=ch1
2808  CALL glk_plframe(id,kax,kay,chxfmt,chyfmt)
2809 * The Y-range from first plot is preserved
2810  CALL glk_range1(id,yl,yu)
2811  ENDIF
2812 * The X-range as in histo
2813  xl = dxl
2814  xu = dxu
2815 *
2816  chrp= chrp1//chrx(incr)
2817  WRITE(m_ltx,'(A)') '%=GLK_Plot2: next plot (line) =========='
2818  WRITE(m_ltx,'(A,I10)')'%====HISTOGRAM ID=',id
2819  WRITE(m_ltx,'(A,A70 )') '% ',title
2820  CALL glk_optout(id,ioplog,iopsla,ioperb,iopsc1,iopsc2)
2821  ker = ioperb-1
2822 * Default line type
2823  IF (iopsla .EQ. 2) THEN
2824  chr='C'
2825  ELSE
2826  chr=' '
2827  ENDIF
2828 * User defined line-type
2829  IF (ch2 .EQ. 'B') chr=' '
2830 *...marks in the midle of the bin
2831  IF (ch2 .EQ. '*') chr='*'
2832 *...marks on the right edge of the bin
2833  IF (ch2 .EQ. 'R') chr='R'
2834 *...marks on the left edge of the bin
2835  IF (ch2 .EQ. 'L') chr='L'
2836  IF (ch2 .EQ. 'C') chr='C'
2837 *...various types of lines
2838  IF (chr .EQ. ' ') THEN
2839 *...contour line used for histogram
2840  CALL glk_plkont(kax,kay,nchx,yl,yu,yy,ker,yer)
2841  ELSE IF(chr .EQ. '*' .OR. chr .EQ. 'R'.OR. chr .EQ. 'L') THEN
2842 *...marks on the right/left/midle of the bin
2843  CALL glk_plmark(kax,kay,nchx,yl,yu,yy,ker,yer,chmark,chr,chrp,chre)
2844  ELSE IF(chr .EQ. 'C') THEN
2845 *...slanted (dotted) line in plotting non-MC functions
2846  CALL glk_plcirc(kax,kay,nchx,yl,yu,yy)
2847  ENDIF
2848 *------------------------------!
2849 * ENDing !
2850 *------------------------------!
2851  WRITE(m_ltx,'(2A)') m_bs,'end{picture} % close entire picture '
2852  IF(abs(m_lint) .EQ. 2) THEN
2853  WRITE(m_ltx,'(A)') '%== GLK_Plot2: end of plot =========='
2854  ELSE
2855  WRITE(m_ltx,'(2A)') m_bs,'end{figure}'
2856  ENDIF
2857  RETURN
2858  900 CALL glk_stop1('+++GLK_Plot2: Nonexistig histo, skipped, id= ',id)
2859  END
2860 
2861  SUBROUTINE glk_plframe(id,kax,kay,chxfmt,chyfmt)
2862 * ************************************************
2863  IMPLICIT NONE
2864  INTEGER id,kax,kay
2865  CHARACTER chxfmt*(*),chyfmt*(*)
2866  include 'GLK.h'
2867  SAVE
2868 *---------------------------------------------------
2869  CHARACTER*80 title
2870  DOUBLE PRECISION dxl,dxu,xl,xu,yl,yu
2871  INTEGER icont,i,nchx
2872  DATA icont/0/
2873 *---------------------------------------------------
2874  icont=icont+1
2875  CALL glk_hinbo1(id,title,nchx,dxl,dxu)
2876  xl = dxl
2877  xu = dxu
2878  CALL glk_range1(id,yl,yu)
2879 *
2880  IF(icont .GT. 1) WRITE(m_ltx,'(2A)') m_bs,'newpage'
2881 *------------------------------!
2882 * Header
2883 *------------------------------!
2884  WRITE(m_ltx,'(A)') ' '
2885  WRITE(m_ltx,'(A)') ' '
2886  WRITE(m_ltx,'(A)')
2887  $'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
2888  WRITE(m_ltx,'(A)')
2889  $'%%%%%%%%%%%%%%%%%%%%%%GLK_PlFrame%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
2890  IF(abs(m_lint) .EQ. 2) THEN
2891  WRITE(m_ltx,'(2A)') m_bs,'noindent'
2892  ELSE
2893  WRITE(m_ltx,'(2A)') m_bs,'begin{figure}[!ht]'
2894  WRITE(m_ltx,'(2A)') m_bs,'centering'
2895  WRITE(m_ltx,'(2A)') m_bs,'htmlimage{scale=1.4}'
2896  ENDIF
2897 *------------------------------!
2898 * General Caption
2899 *------------------------------!
2900  IF(abs(m_lint) .NE. 2) THEN
2901  WRITE(m_ltx,'(6A)')
2902  $ m_bs,'caption{',m_bs,'footnotesize',m_bs,'sf'
2903  DO i=1,m_keytit
2904  WRITE(m_ltx,'(A)') m_titch(i)
2905  ENDDO
2906  WRITE(m_ltx,'(A)') '}'
2907  ENDIF
2908 *------------------------------!
2909 * Frames and labels
2910 *------------------------------!
2911  WRITE(m_ltx,'(A)') '% =========== big frame, title etc. ======='
2912  WRITE(m_ltx,'(4A)') m_bs,'setlength{',m_bs,'unitlength}{0.1mm}'
2913  WRITE(m_ltx,'(2A)') m_bs,'begin{picture}(1600,1500)'
2914  IF( m_lint .LT. 0) THEN
2915 * Big frame usefull for debuging
2916  WRITE(m_ltx,'(4A)')
2917  $ m_bs,'put(0,0){',m_bs,'framebox(1600,1500){ }}'
2918  ENDIF
2919  WRITE(m_ltx,'(A)') '% =========== small frame, labeled axis ==='
2920  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2921  $ m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
2922  WRITE(m_ltx,'(4A,I4,A,I4,A)')
2923  $ m_bs,'put(0,0){',m_bs,'framebox( ',kax,',',kay,'){ }}'
2924  WRITE(m_ltx,'(A)') '% =========== x and y axis ================'
2925  CALL glk_axisx(kax,xl,xu,chxfmt)
2926  CALL glk_axisy(kay,yl,yu,chyfmt)
2927  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}'
2928  $ ,'% end of plotting labeled axis'
2929  END
2930 
2931  SUBROUTINE glk_axisx(kay,yl,yu,chxfmt)
2932 * ***************************************
2933 * plotting x-axis with long and short tips
2934  IMPLICIT NONE
2935  INTEGER kay
2936  DOUBLE PRECISION yl,yu
2937  CHARACTER chxfmt*16
2938  include 'GLK.h'
2939  SAVE
2940 *-------------------------------------------------------
2941  CHARACTER*64 fmt1,fmt2
2942  parameter(fmt1 = '(2A,F8.2,A,F8.2,A,I4,3A)')
2943  parameter(fmt2 = '(2A,F8.2,A,I4,A,F8.2,A,I4,3A)')
2944  DOUBLE PRECISION dy,ddy,ddyl,yy0l,ddys,yy0s,p0s,pds,scmx,p0l,pdl
2945  INTEGER ly,jy,nlt,nts,lex,k,n
2946  DOUBLE PRECISION tipsy(20)
2947 *-------------------------------------------------------
2948  dy= abs(yu-yl)
2949  ly = nint( log10(dy) -0.4999999d0 )
2950  jy = nint(dy/10d0**ly)
2951  ddyl = dy*10d0**(-ly)
2952  IF( jy .EQ. 1) ddyl = 10d0**ly*0.25d0
2953  IF( jy .GE. 2 .AND. jy .LE. 3) ddyl = 10d0**ly*0.5d0
2954  IF( jy .GE. 4 .AND. jy .LE. 6) ddyl = 10d0**ly*1.0d0
2955  IF( jy .GE. 7) ddyl = 10d0**ly*2.0d0
2956  WRITE(m_ltx,'(A)') '% -------GLK_AxisX---- '
2957  WRITE(m_ltx,'(A,I4)') '% JY= ',jy
2958 *-------
2959  nlt = int(dy/ddyl)
2960  nlt = max0(min0(nlt,20),1)+1
2961  yy0l = nint(yl/ddyl+0.5d0)*ddyl
2962  ddys = ddyl/10d0
2963  yy0s = nint(yl/ddys+0.4999999d0)*ddys
2964  p0l = kay*(yy0l-yl)/(yu-yl)
2965  pdl = kay*ddyl/(yu-yl)
2966  p0s = kay*(yy0s-yl)/(yu-yl)
2967  pds = kay*ddys/(yu-yl)
2968  nlt = int(abs(yu-yy0l)/ddyl+0.0000001d0)+1
2969  nts = int(abs(yu-yy0s)/ddys+0.0000001d0)+1
2970  DO n=1,nlt
2971  tipsy(n) =yy0l+ ddyl*(n-1)
2972  ENDDO
2973  WRITE(m_ltx,fmt1)
2974  $ m_bs,'multiput(' ,p0l, ',0)(' ,pdl, ',0){' ,nlt, '}{',
2975  $ m_bs,'line(0,1){25}}',
2976  $ m_bs,'multiput(' ,p0s, ',0)(' ,pds, ',0){' ,nts, '}{',
2977  $ m_bs,'line(0,1){10}}'
2978  WRITE(m_ltx,fmt2)
2979  $ m_bs,'multiput(' ,p0l, ',' ,kay, ')(' ,pdl, ',0){' ,nlt,
2980  $ '}{' ,m_bs, 'line(0,-1){25}}',
2981  $ m_bs,'multiput(' ,p0s, ',' ,kay, ')(' ,pds, ',0){' ,nts,
2982  $ '}{' ,m_bs, 'line(0,-1){10}}'
2983 * ...labeling of axis
2984  scmx = dmax1(dabs(yl),dabs(yu))
2985  lex = nint( log10(scmx) -0.50001)
2986  DO n=1,nlt
2987  k = nint(kay*(tipsy(n)-yl)/(yu-yl))
2988  IF(lex .LE. 3 .AND. lex .GE. -3) THEN
2989 * ...without exponent
2990  WRITE(m_ltx,'(2A,I4,5A,'//chxfmt//',A)')
2991  $ m_bs,'put(',k,',-25){',m_bs,'makebox(0,0)[t]{',
2992  $ m_bs,'Large $ ', tipsy(n), ' $}}'
2993  ELSE
2994 * ...with exponent
2995  WRITE(m_ltx,'(2A,I4,5A,'//chxfmt//',2A,I4,A)')
2996  $ m_bs,'put(' ,k, ',-25){',m_bs,'makebox(0,0)[t]{',
2997  $ m_bs,'Large $ ',
2998  $ tipsy(n)/(10d0**lex),m_bs,'cdot 10^{',lex,'} $}}'
2999  ENDIF
3000  ENDDO
3001  END
3002 
3003  SUBROUTINE glk_axisy(kay,yl,yu,chyfmt)
3004 * ***************************************
3005 * plotting y-axis with long and short tips
3006  IMPLICIT NONE
3007  INTEGER kay
3008  DOUBLE PRECISION yl,yu
3009  CHARACTER chyfmt*16
3010  include 'GLK.h'
3011  SAVE
3012  DOUBLE PRECISION tipsy(20)
3013 *------------------------------------------------------------------
3014  CHARACTER*64 fmt1,fmt2
3015  parameter(fmt1 = '(2A,F8.2,A,F8.2,A,I4,3A)')
3016  parameter(fmt2 = '(2A,I4,A,F8.2,A,F8.2,A,I4,3A)')
3017  INTEGER ly,jy,nlt,nts,lex,n,k
3018  DOUBLE PRECISION ddyl,dy,yy0l,p0l,pdl,pds,scmx,z0l,p0s,yy0s,ddys
3019 *------------------------------------------------------------------
3020  dy= abs(yu-yl)
3021  ly = nint( log10(dy) -0.49999999d0 )
3022  jy = nint(dy/10d0**ly)
3023  ddyl = dy*10d0**(-ly)
3024  IF( jy .EQ. 1) ddyl = 10d0**ly*0.25d0
3025  IF( jy .GE. 2 .AND. jy .LE. 3) ddyl = 10d0**ly*0.5d0
3026  IF( jy .GE. 4 .AND. jy .LE. 6) ddyl = 10d0**ly*1.0d0
3027  IF( jy .GE. 7) ddyl = 10d0**ly*2.0d0
3028  WRITE(m_ltx,'(A)') '% --------GLK_SAxisY------- '
3029  WRITE(m_ltx,'(A,I4)') '% JY= ',jy
3030 *-------
3031  nlt = int(dy/ddyl)
3032  nlt = max0(min0(nlt,20),1)+1
3033  yy0l = nint(yl/ddyl+0.4999999d0)*ddyl
3034  ddys = ddyl/10d0
3035  yy0s = nint(yl/ddys+0.5d0)*ddys
3036  p0l = kay*(yy0l-yl)/(yu-yl)
3037  pdl = kay*ddyl/(yu-yl)
3038  p0s = kay*(yy0s-yl)/(yu-yl)
3039  pds = kay*ddys/(yu-yl)
3040  nlt= int(abs(yu-yy0l)/ddyl+0.0000001d0) +1
3041  nts= int(abs(yu-yy0s)/ddys+0.0000001d0) +1
3042  DO n=1,nlt
3043  tipsy(n) =yy0l+ ddyl*(n-1)
3044  ENDDO
3045 * plotting tics on vertical axis
3046  WRITE(m_ltx,fmt1)
3047  $ m_bs,'multiput(0,' ,p0l, ')(0,' ,pdl ,'){' ,nlt, '}{', m_bs,'line(1,0){25}}',
3048  $ m_bs,'multiput(0,' ,p0s, ')(0,' ,pds, '){' ,nts, '}{', m_bs,'line(1,0){10}}'
3049  WRITE(m_ltx,fmt2)
3050  $ m_bs,'multiput(' ,kay, ',' ,p0l, ')(0,' ,pdl, '){' ,nlt,
3051  $ '}{',m_bs,'line(-1,0){25}}',
3052  $ m_bs,'multiput(' ,kay, ',' ,p0s, ')(0,' ,pds, '){' ,nts,
3053  $ '}{',m_bs,'line(-1,0){10}}'
3054 * ...Zero line if necessary
3055  z0l = kay*(-yl)/(yu-yl)
3056  IF( (z0l .GT. 0d0) .AND. (z0l .LT. float(kay)) )
3057  $ WRITE(m_ltx,'(2A,F8.2,3A,I4,A)') m_bs,'put(0,' ,z0l, '){',m_bs,'line(1,0){' ,kay, '}}'
3058 * ...labeling of axis
3059  scmx = dmax1(dabs(yl),dabs(yu))
3060  lex = nint( log10(scmx) -0.50001d0)
3061  DO n=1,nlt
3062  k = nint(kay*(tipsy(n)-yl)/(yu-yl))
3063  IF(lex .LE. 3 .AND. lex .GE. -3) THEN
3064 * ...without exponent
3065  WRITE(m_ltx,'(2A,I4,5A,'//chyfmt//',A)')
3066  $ m_bs,'put(-25,' ,k, '){',m_bs,'makebox(0,0)[r]{',
3067  $ m_bs,'Large $ ' ,tipsy(n), ' $}}'
3068  ELSE
3069 * ...with exponent
3070  WRITE(m_ltx,'(2A,I4,5A,'//chyfmt//',2A,I4,A)')
3071  $ m_bs,'put(-25,' ,k, '){',m_bs,'makebox(0,0)[r]{',
3072  $ m_bs,'Large $ ',
3073  $ tipsy(n)/(10d0**lex), m_bs,'cdot 10^{' ,lex, '} $}}'
3074  ENDIF
3075  ENDDO
3076  END
3077 
3078  SUBROUTINE glk_plkont(kax,kay,nchx,yl,yu,yy,ker,yer)
3079 */////////////////////////////////////////////////////////////////////////////////////
3080 *// //
3081 *// Plotting contour line for histogram (formely PlHis) //
3082 *// //
3083 */////////////////////////////////////////////////////////////////////////////////////
3084  IMPLICIT NONE
3085  INTEGER kax,kay,nchx,ker
3086  DOUBLE PRECISION yl, yu, yy(*),yer(*),z0l
3087  include 'GLK.h'
3088  SAVE
3089 *---------------------------------------------------
3090  CHARACTER*80 fmt1
3091  INTEGER ix0,iy0,ib,ix1,iy1,ie,ierr,ix2,idy,idx
3092  DOUBLE PRECISION yib
3093 *---------------------------------------------------
3094  WRITE(m_ltx,'(4A,I4,A,I4,A)') m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
3095  WRITE(m_ltx,'(A)') '% ========== plotting primitives =========='
3096 * Color string, optionaly
3097  IF(m_keycol .EQ. 1) THEN
3098  WRITE(m_ltx,'(A)') m_color
3099  m_keycol = 0
3100  ENDIF
3101 *...short macros for vertical/horizontal straight lines
3102  WRITE(m_ltx,'(8A)')
3103  $ m_bs,'newcommand{',m_bs,'x}[3]{',m_bs,'put(#1,#2){', m_bs,'line(1,0){#3}}}'
3104  WRITE(m_ltx,'(8A)')
3105  $ m_bs,'newcommand{',m_bs,'y}[3]{',m_bs,'put(#1,#2){', m_bs,'line(0,1){#3}}}'
3106  WRITE(m_ltx,'(8A)')
3107  $ m_bs,'newcommand{',m_bs,'z}[3]{',m_bs,'put(#1,#2){', m_bs,'line(0,-1){#3}}}'
3108 * error bars
3109  WRITE(m_ltx,'(8A)')
3110  $ m_bs,'newcommand{',m_bs,'e}[3]{', m_bs,'put(#1,#2){',m_bs,'line(0,1){#3}}}'
3111 * Starting point for the line
3112  ix0=0
3113  iy0=0
3114 * Start at Zero line if possible
3115  z0l = kay*(-yl)/(yu-yl)
3116  IF( (z0l .GT. 0d0) .AND. (z0l .LT. float(kay)) ) iy0=z0l
3117  DO ib=1,nchx
3118  yib = yy(ib)
3119  ix1 = nint(kax*(ib-0.00001d0)/nchx) ! new x
3120  iy1 = nint(kay*(yib-yl)/(yu-yl)) ! new y
3121  iy1 = min(max(iy1,-1),kay+1) ! cosmetics
3122  idx = ix1-ix0 ! delta x
3123  idy = iy1-iy0 ! delta y
3124  fmt1 = '(2(2a,i4,a,i4,a,i4,a))'
3125  IF(iy1 .GE. 0 .AND. iy1 .LE. kay) THEN
3126  IF( idy .GE. 0) THEN ! up
3127  WRITE(m_ltx,fmt1) m_bs,'y{',ix0,'}{',iy0,'}{',idy,'}',
3128  $ m_bs,'x{',ix0,'}{',iy1,'}{',idx,'}'
3129  ELSE ! down
3130  WRITE(m_ltx,fmt1) m_bs,'z{',ix0,'}{',iy0,'}{',-idy,'}',
3131  $ m_bs,'x{',ix0,'}{',iy1,'}{',idx,'}'
3132  ENDIF
3133  ENDIF
3134  ix0=ix1
3135  iy0=iy1
3136  IF(ker .EQ. 1) THEN
3137  ix2 = nint(kax*(ib-0.5000d0)/nchx)
3138  ierr = nint(kay*((yy(ib)-yer(ib))-yl)/(yu-yl)) ! bottom of error bar
3139  ie = nint(kay*yer(ib)/(yu-yl)) ! total length of error bar
3140 * Cosmetics
3141  IF(ierr .LT. 0) THEN
3142  ie= ie+ierr
3143  ierr = 0
3144  ENDIF
3145  IF( (ierr+2*ie) .GT. kay) THEN
3146  ie= iabs(kay-ierr)/2
3147  ENDIF
3148  IF( (iy1.GE.0).AND.(iy1.LE. kay).AND.(abs(1d0*ierr).LE.9999d0).AND.(2d0*ie.LE.9999d0) )
3149  $ WRITE(m_ltx,8000) m_bs,ix2,ierr,2*ie
3150  ENDIF
3151  ENDDO
3152 8000 FORMAT(4(a1,2he{,i4,2h}{,i5,2h}{,i4,1h}:1x ))
3153  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}', ' % end of plotting histogram'
3154 * change line-style
3155  m_tline= m_tline+1
3156  IF(m_tline .GT. 2) m_tline=1
3157  END
3158 
3159  SUBROUTINE glk_plmark(kax,kay,nchx,yl,yu,yy,ker,yer,chmark,chr,chr2,chr3)
3160 */////////////////////////////////////////////////////////////////////////////////////
3161 *// //
3162 *// marks in the midle of the bin //
3163 *// //
3164 */////////////////////////////////////////////////////////////////////////////////////
3165  IMPLICIT NONE
3166  INTEGER kax,kay,nchx,ker
3167  DOUBLE PRECISION yl,yu, yy(*),yer(*)
3168  CHARACTER*1 chr
3169  CHARACTER chmark*(*),chr2*(*),chr3*(*)
3170 *---------------------------------------------------
3171  include 'GLK.h'
3172  SAVE
3173  INTEGER ib,ix1,iy1,ierr,ie
3174 *---------------------------------------------------
3175  WRITE(m_ltx,'(4A,I4,A,I4,A)') m_bs,'put(300,250){',m_bs,'begin{picture}( ',kax,',',kay,')'
3176  WRITE(m_ltx,'(A)') '% ===GLK_PlMark: plotting primitives ======'
3177 * Color string, optionaly
3178  IF(m_keycol .EQ. 1) THEN
3179  WRITE(m_ltx,'(A)') m_color
3180  m_keycol = 0
3181  ENDIF
3182 * Plotting symbol
3183  WRITE(m_ltx,'(10A)') m_bs,'newcommand{',m_bs,chr2 , '}[2]{', m_bs,'put(#1,#2){',chmark,'}}'
3184 * Error bar symbol
3185  WRITE(m_ltx,'(10A)')
3186  $ m_bs,'newcommand{',m_bs,chr3 , '}[3]{', m_bs,'put(#1,#2){',m_bs,'line(0,1){#3}}}'
3187 
3188  DO ib=1,nchx
3189  IF(chr .EQ. '*') THEN
3190  ix1 = nint(kax*(ib-0.5000d0)/nchx) ! Midle of bin
3191  ELSEIF(chr .EQ. 'R') THEN
3192  ix1 = nint(kax*(ib*1d0)/nchx) ! Right edge of bin
3193  ELSEIF(chr .EQ. 'L') THEN
3194  ix1 = nint(kax*(ib-1d0)/nchx) ! Left edge of bin
3195  ELSE
3196  WRITE(6,*) '+++++ plamark: wrong line type:',chr
3197  RETURN
3198  ENDIF
3199  iy1 = nint(kay*(yy(ib)-yl)/(yu-yl))
3200  IF(iy1 .GE. 0 .AND. iy1 .LE. kay)
3201  $ WRITE(m_ltx,'(A,A,A,I4,A,I4,A)')
3202  $ m_bs,chr2, '{' ,ix1, '}{' ,iy1, '}'
3203  IF(ker .EQ. 1) THEN
3204  ierr = nint(kay*((yy(ib)-yer(ib))-yl)/(yu-yl)) ! bottom of error bar
3205  ie = nint(kay*yer(ib)/(yu-yl)) ! total length of error bar
3206 * Cosmetics
3207  IF(ierr .LT. 0) THEN
3208  ie= ie+ierr
3209  ierr = 0
3210  ENDIF
3211  IF( (ierr+2*ie) .GT. kay) THEN
3212  ie= iabs(kay-ierr)/2
3213  ENDIF
3214  IF((iy1.GE.0) .AND.(iy1.LE.kay) .AND.(abs(1d0*ierr).LE.9999d0) .AND.(2d0*ie.LE.9999d0))
3215  $ WRITE(m_ltx,'(A,A,A,I4,A,I5,A,I4,A)')
3216  $ m_bs, chr3, '{' ,ix1, '}{' ,ierr, '}{' ,2*ie, '}'
3217  ENDIF
3218  ENDDO
3219  WRITE(m_ltx,'(3A)') m_bs,'end{picture}}',
3220  $ ' % end of plotting histogram'
3221  END
3222 
3223 
3224  SUBROUTINE glk_pltable(Npl,idl,capt,fmt,nch1,incr,npag)
3225 * ******************************************************
3226 * Tables in TeX, up to 9 columns
3227 * Npl = numbers of columns/histograms
3228 * idl(1:Npl) = list of histo id's
3229 * capt(1:Npl+1) = list of captions above each column
3230 * fmt(1:1) = format to print x(i) in first columb,
3231 * h(i) and error he(i) in further columns
3232 * nch1,incr = raws are printet in the sequence
3233 * (h(i),he(i),i=nch1,nbin,incr), nbin is no. of bins.
3234 * npag = 0 no page eject, =1 with page eject
3235 * ******************************************************
3236  IMPLICIT NONE
3237 *--------------- parameters ------------
3238  INTEGER npl,idl(*),nch1,incr,npag
3239  CHARACTER*(*) capt(*)
3240  CHARACTER*(*) fmt(3)
3241 *-------------------------------------------
3242  include 'GLK.h'
3243  SAVE
3244 *---------------------------------------------------
3245  CHARACTER*16 fmt1,fmt2,fmt3
3246  LOGICAL glk_exist
3247  INTEGER i,j,k,n,nchx,nplt,idum,id1,id
3248  INTEGER iopsc1,ioperb,iopsla,iopsc2,ioplog
3249  DOUBLE PRECISION xl,xu,dxl,dxu,xi
3250  DOUBLE PRECISION yyy(m_maxnb),yer(m_maxnb),bi(m_maxnb,9),er(m_maxnb,9)
3251  CHARACTER*80 title
3252  CHARACTER*1 cn(9)
3253  DATA cn /'1','2','3','4','5','6','7','8','9'/
3254 *-----------------------------------------------------------------------------
3255 * Return if histo non-existing or to many columns
3256  IF(.NOT.glk_exist(id)) goto 900
3257  IF(npl .GT. 9 ) goto 901
3258  fmt1 = fmt(1)
3259  fmt2 = fmt(2)
3260  fmt3 = fmt(3)
3261 *
3262 * npack histograms
3263  id1=idl(1)
3264  CALL glk_hinbo1( id1,title,nchx,dxl,dxu)
3265  xl = dxl
3266  xu = dxu
3267  DO n=1,npl
3268  CALL glk_unpak( idl(n),yyy ,' ',idum)
3269  CALL glk_unpak( idl(n),yer ,'ERRO',idum)
3270  DO k=1,nchx
3271  bi(k,n)=yyy(k)
3272  er(k,n)=yer(k)
3273  ENDDO
3274  ENDDO
3275 *------------------------------!
3276 * Header
3277 *------------------------------!
3278  WRITE(m_ltx,'(A)') ' '
3279  WRITE(m_ltx,'(A)') ' '
3280  WRITE(m_ltx,'(A)') '% ========================================='
3281  WRITE(m_ltx,'(A)') '% ============= begin table ==============='
3282  WRITE(m_ltx,'(2A)') m_bs,'begin{table}[!ht]'
3283  WRITE(m_ltx,'(2A)') m_bs,'centering'
3284 *------------------------------!
3285 * Central Caption
3286 *------------------------------!
3287  WRITE(m_ltx,'(4A)') m_bs,'caption{',m_bs,'small'
3288  DO i=1,m_keytit
3289  WRITE(m_ltx,'(A)') m_titch(i)
3290  ENDDO
3291  WRITE(m_ltx,'(A)') '}'
3292 *------------------------------!
3293 * Tabular header
3294 *------------------------------!
3295  WRITE(m_ltx,'(20A)') m_bs,'begin{tabular}
3296  $ {|', ('|c',j=1,npl+1), '||}'
3297 *
3298  WRITE(m_ltx,'(4A)') m_bs,'hline',m_bs,'hline'
3299 *------------------------------!
3300 * Captions in columns
3301 *------------------------------!
3302  WRITE(m_ltx,'(2A)') capt(1),('&',capt(j+1),j=1,npl)
3303 *
3304  WRITE(m_ltx,'(2A)') m_bs,m_bs
3305  WRITE(m_ltx,'(2A)') m_bs,'hline'
3306 *----------------------------------------!
3307 * Table content
3308 * Note that by default RIGHT EDGE of bin is printed, as necessary for
3309 * cumulative distributions, this can be changed with SLAN option
3310 *----------------------------------------!
3311  CALL glk_optout(idl(1),ioplog,iopsla,ioperb,iopsc1,iopsc2)
3312  DO k=nch1,nchx,incr
3313  xi= dxl + (dxu-dxl)*k/(1d0*nchx)
3314  IF(iopsla.eq.2) xi= dxl + (dxu-dxl)*(k-0.5d0)/(1d0*nchx)
3315  IF(ioperb.eq.2) THEN
3316  WRITE(m_ltx,'(A,'//fmt1//','//cn(npl)//'(A,'//fmt2//',A,A,'//fmt3//'), A)')
3317  $ '$', xi, ('$ & $', bi(k,j), m_bs, 'pm', er(k,j), j=1,npl), '$'
3318  WRITE(m_ltx,'(2A)') m_bs,m_bs
3319  ELSE
3320  WRITE(m_ltx,'(A,'//fmt1//','//cn(npl)//'(A,'//fmt2//'), A)')
3321  $ '$', xi, ('$ & $', bi(k,j), j=1,npl), '$'
3322  WRITE(m_ltx,'(2A)') m_bs,m_bs
3323  ENDIF
3324  ENDDO
3325 *------------------------------!
3326 * Ending
3327 *------------------------------!
3328  WRITE(m_ltx,'(4A)') m_bs,'hline',m_bs,'hline'
3329  WRITE(m_ltx,'(2A)') m_bs,'end{tabular}'
3330  WRITE(m_ltx,'(2A)') m_bs,'end{table}'
3331  WRITE(m_ltx,'(A)') '% ============= end table ==============='
3332  WRITE(m_ltx,'(A)') '% ========================================='
3333  IF(npag .NE. 0) WRITE(m_ltx,'(2A)') m_bs,'newpage'
3334 
3335  RETURN
3336  900 CALL glk_retu1('++++ GLK_PlTable: Nonexistig histo id=',id)
3337  RETURN
3338  901 CALL glk_retu1('++++ GLK_PlTable: To many columns Nplt=',nplt)
3339  END
3340 
3341  SUBROUTINE glk_pltable2(Npl,idl,ccapt,mcapt,fmt,chr1,chr2,chr3)
3342 * ***************************************************************
3343 * Tables in TeX, up to 9 columns
3344 * Npl = numbers of columns/histograms
3345 * idl(1:Npl) = list of histo id's
3346 * ccapt(1:Npl+1)= list of column-captions above each column
3347 * mcapt = multicolumn header, none if mcapt=' ',
3348 * fmt(1:1) = format to print x(i) in first columb,
3349 * h(i) and error he(i) in further columns
3350 * chr1 = ' ' normal default, ='S' the Same table continued
3351 * chr2 = ' ' midle of the bin for x(i) in the first column
3352 * = 'R' right edge, ='L' left edge of the bin
3353 * chr3 = ' ' no page eject, ='E' with page eject at the end.
3354 * Furthermore:
3355 * Captions are defined by means of
3356 * CALL GLK_PlCapt(capt) before CALL GLK_PlTable2
3357 * where CHARACTER*80 capt(50) is content of
3358 * caption, line by line, see also comments in GLK_PlCapt routine.
3359 *
3360 * ******************************************************
3361  IMPLICIT NONE
3362 *-------------- parameters--------------
3363  INTEGER npl,idl(*)
3364  CHARACTER*(*) ccapt(*)
3365  CHARACTER*(*) fmt(3)
3366  CHARACTER*1 chr1,chr2,chr3
3367  CHARACTER*(*) mcapt
3368 *----------------------------------------------------------------------
3369  include 'GLK.h'
3370  SAVE
3371 *----------------------------------------------------------------------
3372  CHARACTER*16 fmt1,fmt2,fmt3
3373  LOGICAL glk_exist
3374  INTEGER iopsc1,ioperb,iopsla,iopsc2,ioplog
3375  INTEGER i,j,k,n,idum,id1,id,nchx,nplt
3376  DOUBLE PRECISION xl,xu,xi,dxu,dxl
3377  DOUBLE PRECISION yyy(m_maxnb),yer(m_maxnb),bi(m_maxnb,9),er(m_maxnb,9)
3378  CHARACTER*80 title
3379  CHARACTER*1 cn(9)
3380  INTEGER k1,k2,k3
3381  DATA cn /'1','2','3','4','5','6','7','8','9'/
3382 *----------------------------------------------------------------------
3383 * RETURN if histo non-existing or to many columns
3384  IF(.NOT.glk_exist(id)) goto 900
3385  IF(npl .GT. 9 ) goto 901
3386  fmt1 = fmt(1)
3387  fmt2 = fmt(2)
3388  fmt3 = fmt(3)
3389 *
3390 * unpack histograms
3391  id1 = idl(1)
3392  CALL glk_hinbo1( id1,title,nchx,dxl,dxu)
3393  xl = dxl
3394  xu = dxu
3395  DO n=1,npl
3396  CALL glk_unpak( idl(n),yyy ,' ',idum)
3397  CALL glk_unpak( idl(n),yer ,'ERRO',idum)
3398  DO k=1,nchx
3399  bi(k,n)=yyy(k)
3400  er(k,n)=yer(k)
3401  ENDDO
3402  ENDDO
3403 
3404  IF(chr1 .EQ. ' ' ) THEN
3405 *------------------------------!
3406 * Header
3407 *------------------------------!
3408  WRITE(m_ltx,'(A)') ' '
3409  WRITE(m_ltx,'(A)') ' '
3410  WRITE(m_ltx,'(A)') '% ========================================'
3411  WRITE(m_ltx,'(A)') '% ============ begin table ==============='
3412 *
3413  IF(abs(m_lint) .EQ. 2 ) THEN
3414  WRITE(m_ltx,'(2A)') m_bs,'noindent'
3415  ELSE
3416  WRITE(m_ltx,'(2A)') m_bs,'begin{table}[!ht]'
3417  WRITE(m_ltx,'(2A)') m_bs,'centering'
3418  ENDIF
3419 *------------------------------!
3420 * Central Caption
3421 *------------------------------!
3422  IF(abs(m_lint) .NE. 2 ) THEN
3423  WRITE(m_ltx,'(6A)')
3424  $ m_bs,'caption{',m_bs,'footnotesize',m_bs,'sf'
3425  DO i=1,m_keytit
3426  WRITE(m_ltx,'(A)') m_titch(i)
3427  ENDDO
3428  WRITE(m_ltx,'(A)') '}'
3429  ENDIF
3430 *------------------------------!
3431 * Tabular header
3432 *------------------------------!
3433  WRITE(m_ltx,'(20A)') m_bs,'begin{tabular}
3434  $ {|', ('|c',j=1,npl+1), '||}'
3435  WRITE(m_ltx,'(4A)') m_bs,'hline',m_bs,'hline'
3436 *------------------------------!
3437 * Captions in columns
3438 *------------------------------!
3439  WRITE(m_ltx,'(2A)') ccapt(1),('&',ccapt(j+1),j=1,npl)
3440 *------------------------------!
3441 * Append previous table
3442 *------------------------------!
3443  ELSEIF(chr1 .EQ. 'S' ) THEN
3444  DO i=1,7
3445  backspace(m_ltx)
3446  ENDDO
3447  ELSE
3448  WRITE(*,*) ' ++++ GLK_PlTable2: WRONG chr1 ' ,chr1
3449  ENDIF
3450 
3451  WRITE(m_ltx,'(2A)') m_bs,m_bs
3452  WRITE(m_ltx,'(2A)') m_bs,'hline'
3453 
3454 *------------------------------!
3455 * Optional multicolumn caption
3456 *------------------------------!
3457  IF(mcapt .NE. ' ') THEN
3458  WRITE(m_ltx,'(3A,I2,A)') '& ',m_bs,'multicolumn{',npl,'}{c||}{'
3459  WRITE(m_ltx,'(3A)') ' ',mcapt, ' }'
3460  WRITE(m_ltx,'(2A)') m_bs,m_bs
3461  WRITE(m_ltx,'(2A)') m_bs,'hline'
3462  ENDIF
3463 
3464 *----------------------------------------!
3465 * Table content
3466 * Note that by default RIGHT EDGE of bin is printed, as necessary for
3467 * cumulative distributions, this can be changed with SLAN option
3468 *----------------------------------------!
3469  CALL glk_optout(idl(1),ioplog,iopsla,ioperb,iopsc1,iopsc2)
3470 *
3471 * table printout can be controlled by GLK_SetTabRan(i1,i2,i3)
3472  k1=1
3473  k2=nchx
3474  k3=1
3475  IF( m_keytbr .EQ. 1 ) THEN
3476  k1 = max(k1,m_tabran(1))
3477  k2 = min(k2,m_tabran(2))
3478  k3 = max(k3,m_tabran(3))
3479  m_keytbr = 0
3480  ENDIF
3481 *
3482  DO k=k1,k2,k3
3483  IF(chr2 .EQ. 'R') THEN
3484 * right
3485  xi= dxl + (dxu-dxl)*k/(1d0*nchx)
3486  ELSEIF(chr2 .EQ. 'L') THEN
3487 * left
3488  xi= dxl + (dxu-dxl)*(k-1d0)/(1d0*nchx)
3489  ELSE
3490 * midle
3491  xi= dxl + (dxu-dxl)*(k-0.5d0)/(1d0*nchx)
3492  ENDIF
3493  IF(ioperb.eq.2) THEN
3494  WRITE(m_ltx,'(A,'//fmt1//','//cn(npl)//'(A,'//fmt2//',A,A,'//fmt3//'), A)')
3495  $ '$', xi, ('$ & $', bi(k,j), m_bs, 'pm', er(k,j), j=1,npl), '$'
3496  WRITE(m_ltx,'(2A)') m_bs,m_bs
3497  ELSE
3498  WRITE(m_ltx,'(A,'//fmt1//','//cn(npl)//'(A,'//fmt2//'), A)')
3499  $ '$', xi, ('$ & $', bi(k,j), j=1,npl), '$'
3500  WRITE(m_ltx,'(2A)') m_bs,m_bs
3501  ENDIF
3502  ENDDO
3503 *------------------------------!
3504 * Ending
3505 *------------------------------!
3506  WRITE(m_ltx,'(4A)') m_bs,'hline',m_bs,'hline'
3507  WRITE(m_ltx,'(2A)') m_bs,'end{tabular}'
3508  IF(abs(m_lint) .EQ. 2 ) THEN
3509  WRITE(m_ltx,'(A)') '% ========================================'
3510  ELSE
3511  WRITE(m_ltx,'(2A)') m_bs,'end{table}'
3512  ENDIF
3513  WRITE(m_ltx,'(A)') '% ============= end table =============='
3514  WRITE(m_ltx,'(A)') '% ========================================'
3515  IF(chr3 .EQ. 'E') THEN
3516  WRITE(m_ltx,'(2A)') m_bs,'newpage'
3517  ELSE
3518  WRITE(m_ltx,'(A)') '% ========================================'
3519  ENDIF
3520  RETURN
3521  900 CALL glk_retu1(' ++++ GLK_PlTable2: Nonexistig histo,id= ',id)
3522  RETURN
3523  901 CALL glk_retu1(' ++++ GLK_PlTable2: To many columns Nplt= ',nplt)
3524  END
3525 
3526 
3527  SUBROUTINE glk_wtmon(mode,id,par1,par2,par3)
3528 * ********************************************
3529 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3530 * !!!! It is now replaces by GKL_M package, see below !!!
3531 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3532 *
3533 * Utility program for monitoring M.C. rejection weights.
3534 * ---------------------------------------------------------
3535 * It is backward compatible with WMONIT except:
3536 * (1) for id=-1 one should call as follows:
3537 * GLK_WtMon(-1,id,0d0,1d0,1d0) or skip initialisation completely!
3538 * (2) maximum absolute weight is looked for,
3539 * (3) GLK_Print(-id) prints weight distribution, net profit!
3540 * (4) no restriction id<100 any more!
3541 * ---------------------------------------------------------
3542 * wt is weight, wtmax is maximum weight and rn is random number.
3543 * IF(mode .EQ. -1) then
3544 * initalization if entry id,
3545 * - wtmax is maximum weight used for couting overweighted
3546 * other arguments are ignored
3547 * ELSEIF(mode .EQ. 0) then
3548 * summing up weights etc. for a given event for entry id
3549 * - wt is current weight.
3550 * - wtmax is maximum weight used for couting overweighted
3551 * events with wt>wtmax.
3552 * - rn is random number used in rejection, it is used to
3553 * count no. of accepted (rn < wt/wtmax) and rejected
3554 * (wt > wt/wtmax) events,
3555 * if ro rejection then put rn=0d0.
3556 * ELSEIF(mode .EQ. 1) THEN
3557 * in this mode wmonit repports on accumulated statistics
3558 * - averwt= average weight wt counting all event
3559 * - errela= relative error of averwt
3560 * - nevtot= total number of accounted events
3561 * - nevacc= no. of accepted events (rn < wt/wtmax)
3562 * - nevneg= no. of events with negative weight (wt < 0)
3563 * - nevzer= no. of events with zero weight (wt = 0d0)
3564 * - nevove= no. of overweghted events (wt > wtmax)
3565 * and if you do not want to use cmonit then the value
3566 * the value of averwt is assigned to wt,
3567 * the value of errela is assigned to wtmax and
3568 * the value of wtmax is assigned to rn in this mode.
3569 * ELSEIF(mode .EQ. 2) THEN
3570 * all information defined for entry id defined above
3571 * for mode=2 is just printed of unit nout
3572 * ENDIF
3573 * note that output repport (mode=1,2) is done dynamically just for a
3574 * given entry id only and it may be repeated many times for one id and
3575 * for various id's as well.
3576 * ************************
3577  IMPLICIT NONE
3578  include 'GLK.h'
3579  INTEGER mode,id
3580  DOUBLE PRECISION par1,par2,par3
3581 * locals
3582  INTEGER idg,nevneg,nevzer,nevtot,nevove,nevacc,nbin,lact,ist3,ntot,ist,ist2
3583  DOUBLE PRECISION xl,xu,errela,sswt,averwt,wwmax,swt,wt,wtmax,rn
3584 *---------------------------------------------------------------------------
3585  idg = -id
3586  IF(id .LE. 0) THEN
3587  CALL glk_stop1(' =====> GLK_WtMon: wrong id= ',id)
3588  ENDIF
3589  IF(mode .EQ. -1) THEN
3590 * *******************
3591  nbin = nint(dabs(par3))
3592  IF(nbin .GT. 100) nbin =100
3593  IF(nbin .EQ. 0) nbin =1
3594  xl = par1
3595  xu = par2
3596  IF(xu .LE. xl) THEN
3597  xl = 0d0
3598  xu = 1d0
3599  ENDIF
3600  CALL glk_hadres(idg,lact)
3601  IF(lact .EQ. 0) THEN
3602  CALL glk_book1(idg,' GLK_WtMon $',nbin,xl,xu)
3603  ELSE
3604  WRITE(m_out,*) ' WARNING GLK_WtMon: exists, id= ',id
3605  WRITE( 6,*) ' WARNING GLK_WtMon: exists, id= ',id
3606  ENDIF
3607  ELSEIF(mode .EQ. 0) THEN
3608 * **********************
3609  CALL glk_hadres(idg,lact)
3610  IF(lact .EQ. 0) THEN
3611  WRITE(m_out,*) ' *****> GLK_WtMon: uninitialized, id= ',id
3612  WRITE( 6,*) ' *****> GLK_WtMon: uninitialized, id= ',id
3613  CALL glk_book1(idg,' GLK_WtMon $',1,0d0,1d0)
3614  CALL glk_hadres(idg,lact)
3615  ENDIF
3616  wt =par1
3617  wtmax=par2
3618  rn =par3
3619 * standard entries
3620  CALL glk_fil1(idg,wt,1d0) !!!! <-- principal filling!!!!
3621 * additional goodies
3622  ist = m_index(lact,2)
3623  ist2 = ist+7
3624  ist3 = ist+11
3625 * maximum weight -- maximum by absolute value but keeping sign
3626  m_b(ist3+13) = max( dabs(m_b(ist3+13)) ,dabs(wt))
3627  IF(wt .NE. 0d0) m_b(ist3+13)=m_b(ist3+13) *wt/dabs(wt)
3628 * nevzer,nevove,nevacc
3629  IF(wt .EQ. 0d0) m_b(ist3+10) =m_b(ist3+10) +1d0
3630  IF(wt .GT. wtmax) m_b(ist3+11) =m_b(ist3+11) +1d0
3631  IF(rn*wtmax .LE. wt) m_b(ist3+12) =m_b(ist3+12) +1d0
3632  ELSEIF(mode .GE. 1 .OR. mode .LE. 10) THEN
3633 * *************************************
3634  CALL glk_hadres(idg,lact)
3635  IF(lact .EQ. 0) THEN
3636  CALL glk_stop1(' lack of initialization, id=',id)
3637  ENDIF
3638  ist = m_index(lact,2)
3639  ist2 = ist+7
3640  ist3 = ist+11
3641  ntot = nint(m_b(ist3 +7))
3642  swt = m_b(ist3 +8)
3643  sswt = m_b(ist3 +9)
3644  IF(ntot.LE.0 .OR. swt.EQ.0d0 ) THEN
3645  averwt=0d0
3646  errela=0d0
3647  ELSE
3648  averwt=swt/float(ntot)
3649  errela=sqrt(abs(sswt/swt**2-1d0/float(ntot)))
3650  ENDIF
3651  nevneg = m_b(ist3 +1) !!! it us underflow, xlow=0 assumed!!!
3652  nevzer = m_b(ist3 +10)
3653  nevove = m_b(ist3 +11)
3654  nevacc = m_b(ist3 +12)
3655  wwmax = m_b(ist3 +13)
3656  nevtot = ntot
3657 * Output through parameters
3658  par1 = averwt
3659  par2 = errela
3660  par3 = nevtot
3661  IF(mode .EQ. 2) THEN
3662  par1 = nevacc
3663  par2 = nevneg
3664  par3 = nevove
3665  ELSEIF(mode .EQ. 3) THEN
3666  par1 = nevneg
3667  par2 = wwmax
3668  ENDIF
3669 * no printout for mode <10
3670 * ************************
3671  IF(mode .LE. 9) RETURN
3672  WRITE(m_out,1003) id, averwt, errela, wwmax
3673  WRITE(m_out,1004) nevtot,nevacc,nevneg,nevove,nevzer
3674  IF(mode .LE. 10) RETURN
3675  CALL glk_print(idg)
3676  ELSE
3677 * ****
3678  CALL glk_stop1('+++GLK_WtMon: wrong mode=',mode)
3679  ENDIF
3680 * *****
3681  1003 FORMAT(
3682  $ ' ======================= GLK_WtMon ========================='
3683  $/,' id averwt errela wwmax'
3684  $/, i5, e17.7, f15.9, e17.7)
3685  1004 FORMAT(
3686  $ ' -----------------------------------------------------------'
3687  $/,' nevtot nevacc nevneg nevove nevzer'
3688  $/, 5i12)
3689  END
3690 
3691  SUBROUTINE glk_cumhis(IdGen,id1,id2)
3692 * ************************************
3693 *///////////////////////////////////////////////////////////////////////////
3694 *// Cumulates histogram content starting from UNDERFLOW //
3695 *// and normalizes to the total x-section in NANOBARNS //
3696 *// IdGen is ID of special histogram written by M.C. generator itself //
3697 *// id2. NE. id1 required!!! //
3698 *///////////////////////////////////////////////////////////////////////////
3699 * ***********************************
3700  IMPLICIT NONE
3701  INTEGER idgen,id1,id2
3702 *----------------------------------------------------------------------
3703  include 'GLK.h'
3704  SAVE
3705 *----------------------------------------------------------------------
3706  CHARACTER*80 title
3707  DOUBLE PRECISION x(m_maxnb),er(m_maxnb)
3708  LOGICAL glk_exist
3709  DOUBLE PRECISION swt,sswt,xsec,errel,tmin,tmax
3710  DOUBLE PRECISION xscrnb,erela,wtsup
3711  INTEGER i,nbt,nevt
3712  DOUBLE PRECISION glk_hi,glk_hie
3713 *----------------------------------------------------------------------
3714  IF (glk_exist(id2)) goto 900
3715 *
3716  CALL glk_mgetntot(idgen,nevt)
3717  CALL glk_mgetave( idgen,xscrnb,erela,wtsup)
3718 *
3719  IF(nevt .EQ. 0) goto 901
3720  CALL glk_hinbo1(id1,title,nbt,tmin,tmax)
3721  swt = glk_hi( id1,0) ! UNDERFLOW
3722  sswt = glk_hie(id1,0)**2 ! UNDERFLOW
3723  DO i=1,nbt
3724  swt = swt + glk_hi( id1,i)
3725  sswt = sswt+ glk_hie(id1,i)**2
3726 * note NEVT in error calc. is for the entire sample related
3727 * to the crude x-section XCRU including !!! zero weight events !!!!
3728  xsec = 0d0
3729  errel = 0d0
3730  IF(swt .NE. 0d0 .AND. nevt .NE. 0) THEN
3731  xsec = swt*(xscrnb/nevt)
3732  errel = sqrt(abs(sswt/swt**2-1d0/float(nevt)))
3733  ENDIF
3734  x(i) = xsec
3735  er(i) = xsec*errel
3736  ENDDO
3737 *! store result in id2
3738  CALL glk_book1(id2,title,nbt,tmin,tmax)
3739  CALL glk_pak( id2,x)
3740  CALL glk_pake( id2,er)
3741  CALL glk_idopt(id2,'ERRO')
3742  RETURN
3743  900 WRITE(6,*) '+++++ CUMHIS: ID2 exixsts!!',id2
3744  RETURN
3745  901 WRITE(6,*) '+++++ CUMHIS: EMPTY HISTO ID=',id1
3746  END
3747 
3748 
3749 
3750 
3751  SUBROUTINE glk_renhst(chak,IdGen,id1,id2)
3752 * *****************************************
3753 *///////////////////////////////////////////////////////////////////////////
3754 *// IdGen is ID of special histogram written by M.C. generator itself //
3755 *// This routine RE-NORMALIZES to NANOBARNS or to UNITY //
3756 *// CHAK = 'NB ' normal case [nb] //
3757 *// CHAK = 'NB10' log10 x-scale assumed [nb] //
3758 *// CHAK = 'UNIT' normalization to unity //
3759 *// id2 .NE. id1 required !!! //
3760 *///////////////////////////////////////////////////////////////////////////
3761 * ***********************************
3762  IMPLICIT NONE
3763  CHARACTER*4 chak
3764  INTEGER idgen,id1,id2
3765 *----------------------------------------------------------------------
3766  include 'GLK.h'
3767  SAVE
3768  CHARACTER*80 title
3769  DOUBLE PRECISION xscrnb,erela,wtsup,tmin,tmax
3770  DOUBLE PRECISION swt,fln10,fact
3771  INTEGER i,nbt,nevt
3772  DOUBLE PRECISION glk_hi,glk_hie
3773 *----------------------------------------------------------------------
3774  IF( id2 .eq. id1) goto 900
3775 
3776  CALL glk_mgetntot(idgen,nevt)
3777  CALL glk_mgetave( idgen,xscrnb,erela,wtsup)
3778 *
3779  CALL glk_hinbo1(id1,title,nbt,tmin,tmax)
3780  IF( chak .EQ. 'NB ') THEN
3781  fact = nbt*xscrnb/(nevt*(tmax-tmin))
3782  CALL glk_operat(id1,'+',id1,id2, fact, 0d0)
3783  ELSEIF( chak .EQ. 'NB10') THEN
3784  fln10 = log(10.)
3785  fact = nbt*xscrnb/(nevt*(tmax-tmin)*fln10)
3786  CALL glk_operat(id1,'+',id1,id2, fact, 0d0)
3787  ELSEIF( chak .EQ. 'UNIT') THEN
3788  swt = glk_hi(id1,0)
3789  DO i=1,nbt+1
3790  swt = swt + glk_hi(id1,i)
3791  ENDDO
3792  fact = nbt/((tmax-tmin))/swt
3793  CALL glk_operat(id1,'+',id1,id2, fact, 0d0)
3794  ELSEIF( chak .EQ. 'UN10') THEN
3795  swt = glk_hi(id1,0)
3796  DO i=1,nbt+1
3797  swt = swt + glk_hi(id1,i)
3798  ENDDO
3799  fact = nbt/((tmax-tmin)*log(10.))/swt
3800  CALL glk_operat(id1,'+',id1,id2, fact, 0d0)
3801  ELSEIF( chak .EQ. ' ') THEN
3802  CALL glk_operat(id1,'+',id1,id2, 1d0, 0d0)
3803  ELSE
3804  WRITE(6,*) '+++++ RENHST: wrong chak=',chak
3805  ENDIF
3806 *
3807  RETURN
3808  900 WRITE(6,*) '+++++ RENHST: ID1=ID2=',id1
3809  END
3810 
3811 
3812 *///////////////////////////////////////////////////////////////////////
3813 *// New Weight Motoring ToolBox
3814 *// (replacement for WTmonit etc.)
3815 *//
3816 *// The tool to monitor very precisely the average weigh
3817 *// and other features of the weight distribution.
3818 *// Note that in principle we are vitaly interested in three parts
3819 *// of the weight distribution:
3820 *// Underflow (-infty,0)
3821 *// Regular (0, WTmax)
3822 *// Overflow (WTmax,+infty)
3823 *// with special emphasis on events with exactly zero weight WT=0d0.
3824 *// Nevertheless, we split (0, WTmax) range into several bins
3825 *// in order to be able to visualise the weight distribution.
3826 *// (Using stardard tools for histogram)
3827 *//
3828 *//
3829 *///////////////////////////////////////////////////////////////////////
3830  SUBROUTINE glk_mbook(idm,title,nnchx,WTmax)
3831 * ******************************************
3832 *///////////////////////////////////////////////////////////////////////
3833 *//
3834 *// Booking one entry. Note it is not an ordinary histogram!!!
3835 *// It works just like GLK_Book1 except that it
3836 *// has internaly negative id and x_minimum is always zero.
3837 *//
3838 *///////////////////////////////////////////////////////////////////////
3839  IMPLICIT NONE
3840  include 'GLK.h'
3841  SAVE
3842  INTEGER idm
3843  CHARACTER*80 title
3844  DOUBLE PRECISION wtmax
3845 *
3846  LOGICAL glk_exist
3847  INTEGER j,id,nnchx,nchx,lact,lengt2,ist,ist2,ist3
3848  INTEGER iopsc1, iopsc2, ioperb, ioplog, iopsla
3849  INTEGER iflag1, iflag2
3850  INTEGER ityphi
3851  DOUBLE PRECISION xl,xu,ddx
3852 *-------------------------------------------------
3853  CALL glk_initialize
3854  id = -idm
3855  IF(glk_exist(id)) goto 900
3856  ist=m_length
3857  CALL glk_hadres(0,lact)
3858 * Check if there is a free entry in the m_index
3859  IF(lact .EQ. 0) CALL glk_stop1('GLK_Mbook: no space left,id= ',id)
3860  m_index(lact,1)=id
3861  m_index(lact,2)=m_length
3862  m_index(lact,3)=0
3863 * ---------- limits
3864  CALL glk_copch(title,m_titlc(lact))
3865  nchx =nnchx
3866  IF(nchx .GT. m_maxnb)
3867  $ CALL glk_stop1(' GLK_Mbook: Too many bins ,id= ',id)
3868  xl = 0d0
3869  xu = wtmax
3870 * ---------- title and bin content ----------
3871  lengt2 = m_length +2*nchx +m_buf1+1
3872  IF(lengt2 .GE. m_lenmb)
3873  $ CALL glk_stop1('GLK_Mbook:too litle storage, m_LenmB= ',m_lenmb)
3874 *
3875  DO j=m_length+1,lengt2+1
3876  m_b(j) = 0d0
3877  ENDDO
3878  m_length=lengt2
3879 *... default flags
3880  ioplog = 1
3881  iopsla = 1
3882  ioperb = 1
3883  iopsc1 = 1
3884  iopsc2 = 1
3885  iflag1 =
3886  $ ioplog+10*iopsla+100*ioperb+1000*iopsc1+10000*iopsc2
3887  ityphi = 3 !!!! <-- new type of histo !!!!
3888  iflag2 = ityphi
3889 * examples of decoding flags
3890 * id = nint(m_b(ist+2)-9d0-9d12)/10
3891 * iflag1 = nint(m_b(ist+3)-9d0-9d12)/10
3892 * ioplog = mod(iflag1,10)
3893 * iopsla = mod(iflag1,100)/10
3894 * ioperb = mod(iflag1,1000)/100
3895 * iopsc1 = mod(iflag1,10000)/1000
3896 * iopsc2 = mod(iflag1,100000)/10000
3897 * iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
3898 * ityphi = mod(iflag2,10)
3899 *--------- buffer -----------------
3900 * header
3901  m_b(ist +1) = 9999999999999d0
3902  m_b(ist +2) = 9d12 + id*10 +9d0
3903  m_b(ist +3) = 9d12 + iflag1*10 +9d0
3904  m_b(ist +4) = 9d12 + iflag2*10 +9d0
3905 * dummy vertical scale
3906  m_b(ist +5) = -100d0
3907  m_b(ist +6) = 100d0
3908 * pointer used to speed up search of histogram address
3909  m_b(ist +7) = 0d0
3910 * information on binning
3911  ist2 = ist+7
3912  m_b(ist2 +1) = nchx
3913  m_b(ist2 +2) = xl
3914  m_b(ist2 +3) = xu
3915  ddx = xu-xl
3916  IF(ddx .EQ. 0d0)
3917  $ CALL glk_stop1(' GLK_Mbook: xl=xu, id= ',id)
3918  m_b(ist2 +4) = float(nchx)/ddx
3919 *
3920 * under/over-flow etc.
3921  ist3 = ist+11
3922  DO j=1,13
3923  m_b(ist3 +j)=0d0
3924  ENDDO
3925  RETURN
3926 *----------------
3927  900 CALL glk_retu1(' WARNING GLK_Mbook: already exists id= ', id)
3928  END
3929 
3930 
3931  SUBROUTINE glk_mfill(idm,Wtm,rn)
3932 * ********************************
3933 *///////////////////////////////////////////////////////////////////////
3934 *//
3935 *// filling of M-subpackage entry
3936 *// simillar as fil1 for 1-dim histo but the storage for error
3937 *// is now used to store sum for 'partial averages' <wt-xlowedge>
3938 *//
3939 *///////////////////////////////////////////////////////////////////////
3940  IMPLICIT NONE
3941  INTEGER idm
3942  DOUBLE PRECISION wtm,rn
3943  include 'GLK.h'
3944  SAVE
3945  INTEGER id
3946  INTEGER lact, ist, ist2, ist3, iflag2, nchx, ityphi
3947  INTEGER iposx1,ipose1, kposx1, kpose1, kx
3948  DOUBLE PRECISION wt, deltx, factx, xlowedge
3949  DOUBLE PRECISION xu, xl, x1, wtmax
3950 *---------------------------------
3951  id = -idm
3952  wt = wtm
3953  CALL glk_hadres(id,lact)
3954 * exit for non-existig histo
3955  IF(lact .EQ. 0)
3956  $ CALL glk_stop1('+++GLK_Mfill: nonexisting id= ',id)
3957 
3958  ist = m_index(lact,2)
3959  ist2 = ist+7
3960  ist3 = ist+11
3961 * one-dim. histo only
3962  iflag2 = nint(m_b(ist+4)-9d0-9d12)/10
3963  ityphi = mod(iflag2,10)
3964  IF(ityphi .NE. 3) CALL glk_stop1('+++GLK_Mfill: wrong id= ',id)
3965  x1 = wt
3966  m_index(lact,3)=m_index(lact,3)+1
3967 * for standard average of x=Wt and its error
3968  m_b(ist3 +7) =m_b(ist3 +7) +1
3969  m_b(ist3 +8) =m_b(ist3 +8) +x1
3970  m_b(ist3 +9) =m_b(ist3 +9) +x1*x1
3971 * filling coordinates
3972  nchx = m_b(ist2 +1)
3973  xl = m_b(ist2 +2) !!<--- It was set to zero in book!!!
3974  xu = m_b(ist2 +3)
3975  wtmax = xu
3976  factx = m_b(ist2 +4) ! (fact=nchx/(xu-xl)
3977  deltx = 1d0/factx
3978  IF(x1 .LT. xl) THEN
3979 * (U)nderflow
3980  iposx1 = ist3 +1
3981  ipose1 = ist3 +4
3982  m_b(iposx1) = m_b(iposx1) +1d0
3983  m_b(ipose1) = m_b(ipose1) +wt
3984  ELSEIF(x1 .GT. xu) THEN
3985 * (O)verflow
3986  iposx1 = ist3 +3
3987  ipose1 = ist3 +6
3988  kposx1 = 0
3989  m_b(iposx1) = m_b(iposx1) +1d0
3990  m_b(ipose1) = m_b(ipose1) +(wt- wtmax)
3991  ELSE
3992 * All of (R)egular range (0,WtMax) in one bin
3993  iposx1 = ist3 +2
3994  ipose1 = ist3 +5
3995  m_b(iposx1) = m_b(iposx1) +1d0
3996  m_b(ipose1) = m_b(ipose1) +wt
3997 * (R)egular bin, the ACTUAL one
3998  kx = (x1-xl)*factx+1d0
3999  kx = min( max(kx,1) ,nchx)
4000  kposx1 = ist +m_buf1+kx
4001  kpose1 = ist +m_buf1+nchx+kx
4002  xlowedge = deltx*(kx-1)
4003  m_b(kposx1) = m_b(kposx1) +1d0
4004  m_b(kpose1) = m_b(kpose1) +(wt-xlowedge)
4005  ENDIF
4006 *--------------------------------
4007 * Additional goodies:
4008 * maximum weight -- maximum by absolute value but keeping sign
4009  m_b(ist3+13) = max( dabs(m_b(ist3+13)) ,dabs(wt))
4010  IF(wt .NE. 0d0) m_b(ist3+13)=m_b(ist3+13) *wt/dabs(wt)
4011 * nevzer,nevove,nevacc
4012  IF(wt .EQ. 0d0) m_b(ist3+10) =m_b(ist3+10) +1d0
4013  IF(wt .GT. wtmax) m_b(ist3+11) =m_b(ist3+11) +1d0
4014  IF(rn*wtmax .LE. wt) m_b(ist3+12) =m_b(ist3+12) +1d0
4015 *---
4016  END !GLK_Mfill
4017 
4018 
4019  SUBROUTINE glk_mgetall(idm,
4020  $ avewt,erela, wtsup, avund, avove,
4021  $ ntot,nacc,nneg,nove,nzer)
4022 * ***************************************************************
4023 *///////////////////////////////////////////////////////////////////////
4024 *//
4025 *// Get all statistics out
4026 *//
4027 *///////////////////////////////////////////////////////////////////////
4028  IMPLICIT NONE
4029  INTEGER idm
4030  DOUBLE PRECISION avewt,erela, wtsup, avund, avove
4031  INTEGER ntot,nacc,nneg,nove,nzer
4032  include 'GLK.h'
4033  SAVE
4034  INTEGER id,ist,ist2,ist3,lact
4035  DOUBLE PRECISION swt,sswt
4036 *--------------------
4037  id= -idm
4038  CALL glk_hadres(id,lact)
4039  IF(lact .EQ. 0)
4040  $ CALL glk_stop1('GLK_MgetAll:lack of initialization, id=',id)
4041  ist = m_index(lact,2)
4042  ist2 = ist+7
4043  ist3 = ist+11
4044  ntot = nint(m_b(ist3 +7))
4045  swt = m_b(ist3 +8)
4046  sswt = m_b(ist3 +9)
4047  IF(ntot.LE.0 .OR. swt.EQ.0d0 ) THEN
4048  avewt=0d0
4049  erela=0d0
4050  ELSE
4051  avewt=swt/dfloat(ntot)
4052  erela=sqrt(abs(sswt/swt**2-1d0/float(ntot)))
4053  ENDIF
4054  wtsup = m_b(ist3 +13)
4055  avund = m_b(ist3 +4)/ntot
4056  avove = m_b(ist3 +6)/ntot
4057  nneg = m_b(ist3 +1) ! NB. it is underflow
4058  nzer = m_b(ist3 +10)
4059  nove = m_b(ist3 +11)
4060  nacc = m_b(ist3 +12)
4061 *-----------------------------
4062 * WRITE(m_out,1003) idm, AveWt, ERela, WtSup
4063 * WRITE(m_out,1004) Ntot,Nacc,Nneg,Nove,Nzer
4064 * 1003 FORMAT(
4065 * $ ' ======================= GLK_Mget =========================='
4066 * $/,' id AveWt ERela WtSup'
4067 * $/, i5, e17.7, f15.9, e17.7)
4068 * 1004 FORMAT(
4069 * $ ' -----------------------------------------------------------'
4070 * $/,' Ntot Nacc Nneg Nove Nzer'
4071 * $/, 5i12)
4072 *------------------------------
4073  END
4074 
4075  SUBROUTINE glk_mgetntot(id,Ntot)
4076 *///////////////////////////////////////////////////////////////////////
4077 *//
4078 *// Get Ntotal only
4079 *//
4080 *///////////////////////////////////////////////////////////////////////
4081  IMPLICIT NONE
4082  include 'GLK.h'
4083  SAVE
4084  INTEGER idm,id
4085  DOUBLE PRECISION avewt, erela, wtsup, avund, avove
4086  INTEGER ntot, nacc, nneg, nove, nzer
4087 *--------------------
4088  CALL glk_mgetall(id,
4089  $ avewt,erela, wtsup, avund, avove,
4090  $ ntot,nacc,nneg,nove,nzer)
4091  END
4092 
4093  SUBROUTINE glk_mgetave(id,AveWt,ERela,WtSup)
4094 *///////////////////////////////////////////////////////////////////////
4095 *//
4096 *// Get averages only and highest weight
4097 *//
4098 *///////////////////////////////////////////////////////////////////////
4099  IMPLICIT NONE
4100  include 'GLK.h'
4101  SAVE
4102  INTEGER idm,id
4103  DOUBLE PRECISION avewt, erela, wtsup, avund, avove
4104  INTEGER ntot, nacc, nneg, nove, nzer
4105 *--------------------
4106  CALL glk_mgetall(id,
4107  $ avewt,erela, wtsup, avund, avove,
4108  $ ntot,nacc,nneg,nove,nzer)
4109  END
4110 
4111  SUBROUTINE glk_mprint(idm)
4112 *///////////////////////////////////////////////////////////////////////
4113 *//
4114 *// Printout
4115 *// Note that bin errors have now changed meaning
4116 *//
4117 *///////////////////////////////////////////////////////////////////////
4118  IMPLICIT NONE
4119  include 'GLK.h'
4120  SAVE
4121  INTEGER idm,id
4122  id= -idm
4123  CALL glk_print(id)
4124  END
4125 
4126 *//////////////////////////////////////////////////////////////////////////////
4127 *//////////////////////////////////////////////////////////////////////////////
4128 *//////////////////////////////////////////////////////////////////////////////
4129 *// //
4130 *// Setters and Getters //
4131 *// //
4132 *//////////////////////////////////////////////////////////////////////////////
4133 
4134  SUBROUTINE glk_clone1(id1,id2,title2)
4135 *//////////////////////////////////////////////////////////////////////////////
4136 *// Clone 1-dimensional histo with onl bining and new title //
4137 *//////////////////////////////////////////////////////////////////////////////
4138  CHARACTER*80 title1, title2, title3
4139  INTEGER i,nb
4140  DOUBLE PRECISION xmin,xmax
4141 
4142  CALL glk_hinbo1(id1,title1,nb,xmin,xmax)
4143  CALL glk_copch(title2,title3)
4144  CALL glk_book1(id2,title3,nb,xmin,xmax)
4145 
4146  END
4147 
4148  SUBROUTINE glk_ymimax(id,wmin,wmax)
4149 *//////////////////////////////////////////////////////////////////////////////
4150 *// //
4151 *//////////////////////////////////////////////////////////////////////////////
4152  IMPLICIT NONE
4153  INTEGER id
4154  DOUBLE PRECISION wmin,wmax
4155 
4156  CALL glk_yminim(id,wmin)
4157  CALL glk_ymaxim(id,wmax)
4158  END
4159 
4160 
4161  SUBROUTINE glk_plset(ch,xx)
4162 *//////////////////////////////////////////////////////////////////////////////
4163 *// //
4164 *// Old style seter, sets type of the ploting mark //
4165 *// //
4166 *//////////////////////////////////////////////////////////////////////////////
4167  IMPLICIT NONE
4168  CHARACTER*4 ch
4169  DOUBLE PRECISION xx
4170  include 'GLK.h'
4171  SAVE
4172 *----------------------------------
4173  IF(ch .EQ. 'DMOD') THEN
4174  m_tline = nint(xx)
4175  ENDIF
4176  END
4177 
4178  SUBROUTINE glk_setnout(ilun)
4179 *//////////////////////////////////////////////////////////////////////////////
4180 *// //
4181 *// Sets output unit number //
4182 *// //
4183 *//////////////////////////////////////////////////////////////////////////////
4184  IMPLICIT NONE
4185  include 'GLK.h'
4186  SAVE
4187  INTEGER ilun
4188 
4189  CALL glk_initialize
4190  m_out=ilun
4191  END
4192 
4193  SUBROUTINE glk_getymin(id,ymin)
4194 *//////////////////////////////////////////////////////////////////////////////
4195 *// Sets vertical scale //
4196 *//////////////////////////////////////////////////////////////////////////////
4197  IMPLICIT NONE
4198  include 'GLK.h'
4199  INTEGER id
4200  DOUBLE PRECISION ymin
4201  INTEGER lact,ist
4202 *
4203  CALL glk_hadres(id,lact)
4204  IF(lact .EQ. 0) RETURN
4205  ist= m_index(lact,2)
4206  ymin = m_b(ist+5)
4207  END
4208 
4209  SUBROUTINE glk_getymax(id,ymax)
4210 *//////////////////////////////////////////////////////////////////////////////
4211 *// Sets vertical scale //
4212 *//////////////////////////////////////////////////////////////////////////////
4213  IMPLICIT NONE
4214  include 'GLK.h'
4215  INTEGER id
4216  DOUBLE PRECISION ymax
4217  INTEGER lact,ist
4218 *
4219  CALL glk_hadres(id,lact)
4220  IF(lact .EQ. 0) RETURN
4221  ist= m_index(lact,2)
4222  ymax = m_b(ist+6)
4223  END
4224 
4225  SUBROUTINE glk_setymin(id,ymin)
4226 *//////////////////////////////////////////////////////////////////////////////
4227 *// Sets vertical scale //
4228 *//////////////////////////////////////////////////////////////////////////////
4229  IMPLICIT NONE
4230  include 'GLK.h'
4231  INTEGER id
4232  DOUBLE PRECISION ymin
4233  INTEGER lact,ist
4234 *
4235  CALL glk_hadres(id,lact)
4236  IF(lact .EQ. 0) RETURN
4237  ist= m_index(lact,2)
4238  m_b(ist+5) = ymin
4239  CALL glk_idopt(id,'YMIN')
4240  END
4241 
4242  SUBROUTINE glk_setymax(id,ymax)
4243 *//////////////////////////////////////////////////////////////////////////////
4244 *// Sets vertical scale //
4245 *//////////////////////////////////////////////////////////////////////////////
4246  IMPLICIT NONE
4247  include 'GLK.h'
4248  INTEGER id
4249  DOUBLE PRECISION ymax
4250  INTEGER lact,ist
4251 *
4252  CALL glk_hadres(id,lact)
4253  IF(lact .EQ. 0) RETURN
4254  ist= m_index(lact,2)
4255  m_b(ist+6) = ymax
4256  CALL glk_idopt(id,'YMAX')
4257  END
4258 
4259 
4260  SUBROUTINE glk_getyminymax(id,ymin,ymax)
4261 *//////////////////////////////////////////////////////////////////////////////
4262 *// Sets vertical scale //
4263 *//////////////////////////////////////////////////////////////////////////////
4264  IMPLICIT NONE
4265  include 'GLK.h'
4266  INTEGER id
4267  DOUBLE PRECISION ymin,ymax
4268 *
4269  CALL glk_getymin(id,ymin)
4270  CALL glk_getymax(id,ymax)
4271  END
4272 
4273  SUBROUTINE glk_setyminymax(id,ymin,ymax)
4274 *//////////////////////////////////////////////////////////////////////////////
4275 *// Sets vertical scale //
4276 *//////////////////////////////////////////////////////////////////////////////
4277  IMPLICIT NONE
4278  include 'GLK.h'
4279  INTEGER id
4280  DOUBLE PRECISION ymin,ymax
4281 *
4282  CALL glk_setymin(id,ymin)
4283  CALL glk_setymax(id,ymax)
4284  END
4285 
4286  SUBROUTINE glk_copyymin(id1,id2)
4287 *//////////////////////////////////////////////////////////////////////////////
4288 *// Sets vertical scale //
4289 *//////////////////////////////////////////////////////////////////////////////
4290  IMPLICIT NONE
4291  include 'GLK.h'
4292  INTEGER id1,id2
4293  DOUBLE PRECISION ymin
4294 *
4295  CALL glk_getymin(id1,ymin)
4296  CALL glk_setymin(id2,ymin)
4297  END
4298 
4299  SUBROUTINE glk_copyymax(id1,id2)
4300 *//////////////////////////////////////////////////////////////////////////////
4301 *// Sets vertical scale //
4302 *//////////////////////////////////////////////////////////////////////////////
4303  IMPLICIT NONE
4304  include 'GLK.h'
4305  INTEGER id1,id2
4306  DOUBLE PRECISION ymax
4307 *
4308  CALL glk_getymax(id1,ymax)
4309  CALL glk_setymax(id2,ymax)
4310  END
4311 
4312  SUBROUTINE glk_setcolor(Color)
4313 *//////////////////////////////////////////////////////////////////////////////
4314 *// //
4315 *// Sets output unit number //
4316 *// //
4317 *//////////////////////////////////////////////////////////////////////////////
4318  IMPLICIT NONE
4319  include 'GLK.h'
4320  CHARACTER*(*) color
4321 *
4322  CALL glk_copch(color,m_color)
4323 *
4324  m_keycol = 1 !flag set up
4325  END
4326 
4327  SUBROUTINE glk_settabran(i1,i2,i3)
4328 *//////////////////////////////////////////////////////////////////////////////
4329 *// //
4330 *// Sets table range for taple printout in GKL_PlTable2 //
4331 *// i1,i2,i3 are lower limit, upper limit and increment //
4332 *// //
4333 *//////////////////////////////////////////////////////////////////////////////
4334  IMPLICIT NONE
4335  include 'GLK.h'
4336  INTEGER i1,i2,i3
4337 *
4338  m_keytbr = 1 !flag set up
4339  m_tabran(1) = i1
4340  m_tabran(2) = i2
4341  m_tabran(3) = i3
4342  END
4343 
4344  SUBROUTINE glk_getnb(id,Nb)
4345 *//////////////////////////////////////////////////////////////////////////////
4346 *// //
4347 *//////////////////////////////////////////////////////////////////////////////
4348  IMPLICIT NONE
4349  include 'GLK.h'
4350  INTEGER id,nb
4351 * local
4352  CHARACTER*80 title
4353  INTEGER lact,ist,ist2
4354 *-------------
4355  CALL glk_hadres(id,lact)
4356  IF(lact .EQ. 0) THEN
4357  CALL glk_stop1('+++STOP in GLK_GetNb: wrong id=',id)
4358  ENDIF
4359  ist = m_index(lact,2)
4360  ist2 = ist+7
4361  nb = m_b(ist2 +1)
4362  END
4363 
4364  SUBROUTINE glk_getbin(id,ib,Bin)
4365 *//////////////////////////////////////////////////////////////////////////////
4366 *// //
4367 *// getting out bin content //
4368 *// //
4369 *//////////////////////////////////////////////////////////////////////////////
4370  IMPLICIT NONE
4371  include 'GLK.h'
4372  INTEGER id,ib
4373  DOUBLE PRECISION bin,glk_hi
4374 
4375  bin = glk_hi(id,ib)
4376  END
4377 
4378 
4379 *//////////////////////////////////////////////////////////////////////////////
4380 *// //
4381 *// End of CLASS GLK //
4382 *//////////////////////////////////////////////////////////////////////////////
4383