C++ Interface to Tauola
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
24708000 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
25417000 FORMAT(4(a1,2hr{,i4,2h}{,i4,1h}:1x ))
25428000 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)
26377000 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
31528000 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