SND@LHC Software
Loading...
Searching...
No Matches
mphistab.f90
Go to the documentation of this file.
1
2! Code converted using TO_F90 by Alan Miller
3! Date: 2012-03-16 Time: 11:07:45
4
89
90! *************************** Histograms ******************************
91
92SUBROUTINE hmpdef(ih,xa,xb,text) ! book, reset histogram
93 USE mpdef
94
95 IMPLICIT NONE
96 INTEGER(mpi) :: i
97 INTEGER(mpi) :: iha
98 INTEGER(mpi) :: ihb
99 INTEGER(mpi) :: ihc
100 INTEGER(mpi) :: ix
101 INTEGER(mpi) :: j
102 INTEGER(mpi) :: lun
103 INTEGER(mpi) :: lunw
104 INTEGER(mpi) :: nbin
105 INTEGER(mpi) :: nn
106 REAL(mps) :: x
107 REAL(mps) :: xcent
108 REAL(mps) :: xmean
109 REAL(mps) :: xsigm
110 ! book millepede histogram, 120 bins
111
112 INTEGER(mpi), INTENT(IN) :: ih
113 REAL(mps), INTENT(IN) :: xa
114 REAL(mps), INTENT(IN) :: xb
115 CHARACTER (LEN=*), INTENT(IN) :: text
116 INTEGER(mpi), PARAMETER :: numhis=15
117 INTEGER(mpi) :: inhist(120,numhis)
118 INTEGER(mpi) ::jnhist(5,numhis)
119 INTEGER(mpi) ::khist(numhis)
120 REAL(mps) :: fnhist(120,numhis)
121 equivalence(inhist(1,1),fnhist(1,1))
122 INTEGER(mpi) :: kvers(numhis)
123 REAL(mps) :: xl(6,numhis)
124 REAL(mpd):: dl(2,numhis)
125 CHARACTER (LEN=60):: htext(numhis)
126 SAVE
127 DATA khist/numhis*0/,lun/7/
128 ! ...
129 IF(ih <= 0.OR.ih > numhis) RETURN
130 ! IF(XA.EQ.XB) RETURN
131 DO i=1,120
132 inhist(i,ih)=0
133 END DO
134 DO j=1,5
135 jnhist(j,ih)=0
136 END DO
137 xl(1,ih)=xa
138 xl(2,ih)=xb
139 xl(3,ih)=0.0
140 IF(xa /= xb) xl(3,ih)=120.0/(xb-xa)
141 xl(6,ih)=0.5*(xa+xb) ! center
142 IF(khist(ih) == 0) THEN
143 kvers(ih)=0
144 ELSE
145 kvers(ih)=kvers(ih)+1
146 END IF
147 khist(ih)=1 ! flt.pt. (lin)
148 htext(ih)=text
149 dl(1,ih)=0.0_mpd
150 dl(2,ih)=0.0_mpd
151 RETURN
152
153 entry hmpldf(ih,text) ! book, reset log histogram
154 IF(ih <= 0.OR.ih > numhis) RETURN
155 DO i=1,120
156 inhist(i,ih)=0
157 END DO
158 DO j=1,5
159 jnhist(j,ih)=0
160 END DO
161 IF(khist(ih) == 0) THEN
162 kvers(ih)=0
163 ELSE
164 kvers(ih)=kvers(ih)+1
165 END IF
166 khist(ih)=2 ! integer log
167 htext(ih)=text
168 xl(1,ih)=0.0
169 xl(2,ih)=6.0
170 RETURN
171
172 entry hmpent(ih,x) ! entry flt.pt.
173 IF(ih <= 0.OR.ih > numhis) RETURN
174 IF(khist(ih) /= 1) RETURN
175 IF(jnhist(4,ih) >= 2147483647) RETURN
176 jnhist(4,ih)=jnhist(4,ih)+1 ! count
177 IF(jnhist(4,ih) <= 120) THEN
178 fnhist(jnhist(4,ih),ih)=x ! store value
179 IF(jnhist(4,ih) == 120) THEN
180 CALL hmpmak(inhist(1,ih),fnhist(1,ih),jnhist(1,ih), xl(1,ih),dl(1,ih))
181 END IF
182 RETURN
183 END IF
184 ! IF(JNHIST(1,IH)+JNHIST(2,IH)+JNHIST(3,IH).EQ.0) THEN
185 ! XL(4,IH)=X
186 ! XL(5,IH)=X
187 ! END IF
188 i=int(1.0+xl(3,ih)*(x-xl(1,ih)),mpi) ! X - Xmin
189 j=2
190 IF(i < 1) j=1
191 IF(i > 120) j=3
192 jnhist(j,ih)=jnhist(j,ih)+1
193 xl(4,ih)=min(xl(4,ih),x)
194 xl(5,ih)=max(xl(5,ih),x)
195 IF(j /= 2) RETURN
196 inhist(i,ih)=inhist(i,ih)+1
197 dl(1,ih)=dl(1,ih)+ x-xl(6,ih)
198 dl(2,ih)=dl(2,ih)+(x-xl(6,ih))**2
199 RETURN
200
201 entry hmplnt(ih,ix) ! entry integer
202 IF(ih <= 0.OR.ih > numhis) RETURN
203 IF(khist(ih) /= 2) RETURN
204 IF(jnhist(1,ih) >= 2147483647) RETURN
205 IF(ix <= 0) THEN
206 jnhist(1,ih)=jnhist(1,ih)+1
207 ELSE
208 IF(jnhist(4,ih) == 0) jnhist(4,ih)=ix
209 IF(jnhist(5,ih) == 0) jnhist(5,ih)=ix
210 jnhist(4,ih)=min(jnhist(4,ih),ix)
211 jnhist(5,ih)=max(jnhist(5,ih),ix)
212 i=int(1.0+20.0*log10(real(ix,mps)),mpi)
213 j=2
214 IF(i < 1) j=1
215 IF(i > 120) j=3
216 IF(j == 2) inhist(i,ih)=inhist(i,ih)+1
217 jnhist(j,ih)=jnhist(j,ih)+1
218 END IF
219 RETURN
220
221 entry hmprnt(ih) ! print, content vert
222 IF(ih == 0) THEN
223 iha=1
224 ihb=numhis
225 ELSE
226 IF(ih <= 0.OR.ih > numhis) RETURN
227 iha=ih
228 ihb=ih
229 END IF
230 DO ihc=iha,ihb
231 IF(khist(ihc) /= 0) THEN
232 IF(khist(ihc) == 1) THEN
233 CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), &
234 xl(1,ihc),dl(1,ihc))
235 END IF
236 nn=jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc)
237 IF(nn /= 0.OR.khist(ihc) == 3) THEN
238 WRITE(*,111)
239111 FORMAT(' ______',2('______________________________'))
240 IF(kvers(ihc) == 1) THEN
241 WRITE(*,*) 'Histogram',ihc,': ',htext(ihc)
242 ELSE
243 WRITE(*,*) 'Histogram',ihc,'/',kvers(ihc),': ',htext(ihc)
244 END IF
245 IF(khist(ihc) == 1) THEN
246 WRITE(*,*) ' Out_low inside out_high = ', (jnhist(j,ihc),j=1,3)
247 ELSE IF(khist(ihc) == 2) THEN
248 WRITE(*,*) ' 0_or_negative inside above_10^6 = ', &
249 (jnhist(j,ihc),j=1,3)
250 END IF
251 IF(khist(ihc) == 3) THEN
252 CALL pfvert(120,fnhist(1,ihc))
253 END IF
254 IF(jnhist(2,ihc) /= 0) THEN ! integer content
255 CALL pivert(120,inhist(1,ihc))
256 IF(khist(ihc) == 1) THEN
257 CALL psvert(xl(1,ihc),xl(2,ihc))
258 ELSE IF(khist(ihc) == 2) THEN
259 CALL psvert(0.0,6.0)
260 END IF
261 END IF
262 IF(khist(ihc) == 1) THEN
263 WRITE(*,*) ' Min and Max are',xl(4,ihc),xl(5,ihc)
264 IF(jnhist(2,ihc) > 1) THEN
265 xmean=real(xl(6,ihc)+dl(1,ihc)/real(jnhist(2,ihc),mps),mps)
266 xcent=0.5*(xl(1,ihc)+xl(2,ihc))
267 xsigm=real((dl(2,ihc)-dl(1,ihc)**2/real(jnhist(2,ihc),mps)),mps)
268 xsigm=sqrt(xsigm/real(jnhist(2,ihc)-1,mps))
269 WRITE(*,*) ' Mean and sigma are', xmean,' +-',xsigm
270 END IF
271 ELSE IF(khist(ihc) == 2) THEN
272 WRITE(*,*) ' Plot of log10 of entries. Min and Max are', &
273 jnhist(4,ihc),jnhist(5,ihc)
274 END IF
275 END IF
276 END IF
277 END DO
278 RETURN
279
280 entry hmplun(lunw) ! unit for output
281 lun=lunw
282 RETURN
283
284 entry hmpwrt(ih) ! write histogram text file
285 IF(lun <= 0) RETURN
286 IF(ih == 0) THEN
287 iha=1
288 ihb=numhis
289 ELSE
290 IF(ih <= 0.OR.ih > numhis) RETURN
291 iha=ih
292 ihb=ih
293 END IF
294
295 DO ihc=iha,ihb ! histogram loop
296 IF(khist(ihc) /= 0) THEN
297 IF(khist(ihc) == 1) THEN
298 CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), &
299 xl(1,ihc),dl(1,ihc))
300 END IF
301 nbin=120
302 WRITE(lun,204) ' '
303 WRITE(lun,201) ihc,kvers(ihc),khist(ihc)
304 WRITE(lun,204) htext(ihc)
305 IF (jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc) == 0 &
306 .AND.xl(1,ihc) == xl(2,ihc)) THEN
307 ! hist is empty and hist range makes no sense
308 ! - cause: hist with 'variable edges' was never filled
309 ! - workaround: make lower and upper edge of hist differ in output
310 WRITE(lun,202) nbin,xl(1,ihc)-0.001,xl(2,ihc)+0.001
311 ELSE
312 WRITE(lun,202) nbin,xl(1,ihc),xl(2,ihc)
313 END IF
314 WRITE(lun,203) (jnhist(j,ihc),j=1,3)
315 WRITE(lun,204) 'bincontent'
316 IF(khist(ihc) == 1.OR.khist(ihc) == 2) THEN
317 CALL kprint(lun,inhist(1,ihc),nbin)
318 ELSE
319 WRITE(lun,219) (fnhist(i,ihc),i=1,nbin)
320 END IF
321
322 IF(khist(ihc) == 1) THEN
323 WRITE(lun,205) xl(4,ihc),xl(5,ihc)
324 ELSE IF(khist(ihc) == 2) THEN
325 WRITE(lun,205) real(jnhist(4,ihc),mps),real(jnhist(5,ihc),mps)
326 END IF
327 IF(khist(ihc) == 1) THEN
328 IF(jnhist(2,ihc) > 1) THEN
329 xmean=real(xl(6,ihc)+dl(1,ihc)/real(jnhist(2,ihc),mps),mps)
330 xcent=0.5*(xl(1,ihc)+xl(2,ihc))
331 xsigm=real((dl(2,ihc)-dl(1,ihc)**2/real(jnhist(2,ihc),mps)),mps)
332 xsigm=sqrt(xsigm/real(jnhist(2,ihc)-1,mps))
333 WRITE(lun,206) xmean,xsigm
334 END IF
335 END IF
336 WRITE(lun,204) 'end of histogram'
337 END IF
338 END DO
339
340201 FORMAT('Histogram ',i4,10x,'version ',i4,10x,'type',i2)
341202 FORMAT(10x,' bins, limits ',i4,2g15.5)
342203 FORMAT(10x,'out-low inside out-high ',3i10)
343204 FORMAT(a)
344205 FORMAT('minmax',2e15.7)
345206 FORMAT('meansigma',2e15.7)
346
347219 FORMAT(4e15.7)
348END SUBROUTINE hmpdef
349
350SUBROUTINE hmpmak(inhist,fnhist,jnhist,xl,dl) ! hist scale from data
351 USE mpdef
352
353 IMPLICIT NONE
354 INTEGER(mpi) :: i
355 INTEGER(mpi) :: j
356 INTEGER(mpi) :: k
357 INTEGER(mpi) :: nn
358 REAL(mps) :: x
359 REAL(mps) :: xa
360 REAL(mps) :: xb
361
362 INTEGER(mpi), INTENT(OUT) :: inhist(120)
363 REAL(mps), INTENT(IN) :: fnhist(120)
364 INTEGER(mpi), INTENT(IN OUT) :: jnhist(5)
365 REAL(mps), INTENT(IN OUT) :: xl(6)
366 REAL(mpd), INTENT(OUT) :: dl(2)
367 REAL(mps) :: cphist(120)
368
369
370
371 SAVE
372 ! ...
373 nn=jnhist(4)
374 ! WRITE(*,*) 'HMPMAK: NN,JNHIST(5)',NN,JNHIST(5)
375 IF(nn == 0.OR.jnhist(5) /= 0) RETURN
376 jnhist(5)=1
377 DO i=1,nn
378 ! WRITE(*,*) 'copy ',I,FNHIST(I)
379 cphist(i)=fnhist(i)
380 END DO
381 CALL heapf(cphist,nn)
382 IF(xl(3) == 0.0) THEN
383 CALL bintab(cphist,nn,xa,xb)
384 xl(1)=xa
385 xl(2)=xb
386 xl(3)=0.0
387 IF(xa /= xb) xl(3)=120.0/(xb-xa)
388 xl(6)=0.5*(xa+xb) ! center
389 END IF
390 xl(4)=cphist( 1)
391 xl(5)=cphist(nn)
392 ! WRITE(*,*) 'XL ',XL
393 DO i=1,nn
394 inhist(i)=0
395 END DO
396 DO k=1,nn
397 x=cphist(k)
398 i=int(1.0+xl(3)*(x-xl(1)),mpi) ! X - Xmin
399 ! WRITE(*,*) 'K,I,X ',K,I,X
400 j=2
401 IF(i < 1) j=1
402 IF(i > 120) j=3
403 jnhist(j)=jnhist(j)+1
404 IF(j == 2) THEN
405 inhist(i)=inhist(i)+1
406 dl(1)=dl(1)+ x-xl(6)
407 dl(2)=dl(2)+(x-xl(6))**2
408 END IF
409 END DO
410END SUBROUTINE hmpmak
411
412
413
414
415SUBROUTINE bintab(tab,n,xa,xb) ! hist scale from data
416 USE mpdef
417
418 IMPLICIT NONE
419 REAL(mps) :: dd
420 REAL(mps) :: dx
421 INTEGER(mpi) :: i
422 INTEGER(mpi) :: iexp
423 INTEGER(mpi) :: ii
424 INTEGER(mpi) :: j
425 INTEGER(mpi) :: m1
426 INTEGER(mpi) :: m2
427 INTEGER(mpi) :: n1
428 INTEGER(mpi) :: n2
429 REAL(mps) :: rat
430 REAL(mps) :: x1
431 REAL(mps) :: x2
432 REAL(mps) :: xx
433 ! Bin limits XA and XB from TAB(N)
434
435 REAL(mps), INTENT(IN) :: tab(n)
436 INTEGER(mpi), INTENT(IN) :: n
437 REAL(mps), INTENT(OUT) :: xa
438 REAL(mps), INTENT(OUT) :: xb
439
440 REAL(mps) :: bin(10)
441 DATA bin/1.0,1.5,2.0,3.0,4.0,5.0,8.0,10.0,15.0,20.0/
442 SAVE
443 ! ...
444
445 CALL heapf(tab,n) ! reduced statistic
446 ! WRITE(*,*) ' '
447 ! WRITE(*,*) 'Sorted ',(TAB(I),I=1,N)
448 IF(n < 100) THEN
449 x1=tab(1)
450 x2=tab(n)
451 ! WRITE(*,*) 'reduced statistic X1 X2 ',X1,X2
452 ELSE ! large statistic
453 m1=int(1.0+0.05*real(n),mpi)
454 m2=int(1.0+0.16*real(n),mpi)
455 x1=tab(m1)-4.0*(tab(m2)-tab(m1))
456 IF(x1 < 0.0.AND.tab(1) >= 0.0) x1=tab(1)
457 x2=tab(n+1-m1)+4.0*(tab(n+1-m1)-tab(n+1-m2))
458 IF(x2 > 0.0.AND.tab(n) <= 0.0) x2=tab(n)
459 ! WRITE(*,*) 'large statistic ',X1,X2
460 ! WRITE(*,*) 'min und max ',TAB(1),TAB(N)
461 IF(x1*tab(1) <= 0.0) x1=0.0
462 IF(x2*tab(n) <= 0.0) x2=0.0
463 ! WRITE(*,*) 'large statistic zero ',X1,X2
464 IF(x1*x2 < 0.0.AND.min(-x1,x2) > 0.6*max(-x1,x2)) THEN
465 xx=max(-x1,x2) ! symmetry
466 x1=-xx
467 x2=+xx
468 ELSE IF(x1*x2 > 0.0.AND. & ! include zero ?
469 abs(min(x1,x2)) < 0.4*abs(max(x1,x2))) THEN
470 IF(x1 < 0.0) THEN
471 x2=0.0
472 ELSE
473 x1=0.0
474 END IF
475 END IF
476 ! WRITE(*,*) 'large statistic ',X1,X2
477 END IF
478 IF(x1 == x2) THEN
479 x1=x1-1.0
480 x2=x2+1.0
481 END IF
482 dx=x2-x1
483 ! WRITE(*,*) 'X1,X2,DX ',X1,X2,DX
484 rat=0.0
485 ii=1
486 DO j=1,11
487 i=j
488 IF(j == 11) i=ii
489 iexp=int(101.0+log10(dx)-log10(6.0*bin(i)),mpi)
490 iexp=iexp-100
491 dd=bin(i)*10.0**iexp
492
493 n1=int(abs(x1)/dd,mpi)
494 IF(x1 < 0.0) n1=-n1
495 IF(real(n1,mps)*dd > x1) n1=n1-1
496 ! WRITE(*,*) 'Bin ',I,N1,N1*DD,X1
497
498 n2=int(abs(x2)/dd,mpi)
499 IF(x2 < 0.0) n2=-n2
500 IF(real(n2,mps)*dd < x2) n2=n2+1
501 ! WRITE(*,*) 'Bin ',I,N2,N2*DD,X2
50210 IF(n2-n1 < 6) THEN
503 IF(n1 /= 0) n1=n1-1
504 IF(n2-n1 < 6.AND.n2 /= 0) n2=n2+1
505 GO TO 10
506 END IF
507 ! WRITE(*,*) 'corrected N1 N2 ',N1,N2
508 xa=sign(real(n1,mps)*dd,x1)
509 xb=sign(real(n2,mps)*dd,x2)
510 ! WRITE(*,*) J,' resulting limits XA XB ',XA,XB
511 IF((x2-x1)/(xb-xa) > rat) THEN
512 ii=i
513 rat=(x2-x1)/(xb-xa)
514 END IF
515 END DO
516! WRITE(*,*) J,' resulting limits XA XB ',XA,XB
517END SUBROUTINE bintab
518
519SUBROUTINE kprint(lun,list,n) ! print integer array
520 USE mpdef
521
522 IMPLICIT NONE
523 INTEGER(mpi) :: i
524 INTEGER(mpi) :: ia
525 INTEGER(mpi) :: ib
526 INTEGER(mpi) :: k
527 INTEGER(mpi) :: ln
528 INTEGER(mpi) :: lp
529 INTEGER(mpi) :: np
530 ! print integer array LIST(N)
531
532 INTEGER(mpi), INTENT(IN OUT) :: lun
533 INTEGER(mpi), INTENT(IN) :: list(n)
534 INTEGER(mpi), INTENT(IN) :: n
535 INTEGER(mpi) :: li(7)
536 DATA li/2,3,4,6,8,9,12/ ! number of characters
537 SAVE
538 ! ...
539 ib=0
54010 ia=ib+1
541 IF(ia > n) RETURN
542 DO k=1,7
543 np=72/li(k)
544 ib=min(ia-1+np,n)
545 IF(k <= 6) THEN
546 lp=10**(li(k)-1)-1 ! maximum positive
547 ln=-lp/10 ! minimum negative
548 DO i=ia,ib
549 IF(list(i) > lp.OR.list(i) < ln) GO TO 20
550 END DO
551 END IF
552 IF(k == 1) THEN
553 WRITE(lun,101) (list(i),i=ia,ib)
554 ELSE IF(k == 2) THEN
555 WRITE(lun,102) (list(i),i=ia,ib)
556 ELSE IF(k == 3) THEN
557 WRITE(lun,103) (list(i),i=ia,ib)
558 ELSE IF(k == 4) THEN
559 WRITE(lun,104) (list(i),i=ia,ib)
560 ELSE IF(k == 5) THEN
561 WRITE(lun,105) (list(i),i=ia,ib)
562 ELSE IF(k == 6) THEN
563 WRITE(lun,106) (list(i),i=ia,ib)
564 ELSE IF(k == 7) THEN
565 WRITE(lun,107) (list(i),i=ia,ib)
566 END IF
567 GO TO 10
56820 CONTINUE
569 END DO
570101 FORMAT(36i2)
571102 FORMAT(24i3)
572103 FORMAT(18i4)
573104 FORMAT(12i6)
574105 FORMAT( 9i8)
575106 FORMAT( 8i9)
576107 FORMAT( 6i12)
577END SUBROUTINE kprint
578
579! ***************************** XY data ****************************
580
581SUBROUTINE gmpdef(ig,ityp,text) ! book, reset XY storage
582 USE mpdef
583
584 IMPLICIT NONE
585 REAL(mps) :: dx
586 REAL(mps) :: dy
587 INTEGER(mpi) :: i
588 INTEGER(mpi) :: iga
589 INTEGER(mpi) :: igb
590 INTEGER(mpi) :: igc
591 INTEGER(mpi) :: j
592 INTEGER(mpi) :: lun
593 INTEGER(mpi) :: lunw
594 INTEGER(mpi) :: n
595 INTEGER(mpi) :: na
596 REAL(mps) :: wght
597 REAL(mps) :: x
598 REAL(mps) :: y
599 REAL(mps) :: y1
600 ! ITYP = 1 X,Y as dots
601 ! = 2 X,Y as line
602 ! = 3 X,Y as line and dots
603 ! = 4 X,Y, DX,DY symbols
604
605 INTEGER(mpi), INTENT(IN) :: ig
606 INTEGER(mpi), INTENT(IN) :: ityp
607 CHARACTER (LEN=*), INTENT(IN) :: text
608 INTEGER(mpi), PARAMETER :: narr=1000
609 REAL(mps) :: array(2,narr)
610 REAL(mps) ::array4(4,narr/2)
611 REAL(mps) ::array1(narr+narr)
612 REAL(mps) ::four(4)
613 equivalence(array(1,1),array4(1,1),array1(1))
614 INTEGER(mpi), PARAMETER :: numgxy=10
615 INTEGER(mpi), PARAMETER :: nlimit=500
616 INTEGER(mpi) :: nstr(numgxy)
617 INTEGER(mpi) ::igtp(numgxy)
618 INTEGER(mpi) ::lvers(numgxy)
619 INTEGER(mpi) ::nst(3,numgxy)
620 REAL(mps) :: xyplws(10,numgxy)
621 INTEGER(mpi) :: jflc(5,numgxy)
622 INTEGER(mpi) ::kflc(5,numgxy)
623 ! JFLC(1,.) = first used index
624 ! JFLC(2,.) = last used index
625 ! JFLC(3,.) = counter of used places
626 ! JFLC(4,.) = counter of ignored
627 ! JFLC(5,.) = limit for JFLC(3)
628 CHARACTER (LEN=60):: gtext(numgxy)
629
630 LOGICAL:: start
631 SAVE
632 DATA start/.true./,lun/7/
633 DATA nstr/numgxy*0/ ! by GF
634 ! ...
635 IF(start) THEN
636 start=.false.
637 CALL stmars ! initialize storage
638 DO i=1,numgxy
639 DO j=1,5
640 jflc(j,i)=0
641 kflc(j,i)=0
642 END DO
643 END DO
644 END IF
645
646 IF(ig < 1.OR.ig > numgxy) RETURN
647 IF(ityp < 1.OR.ityp > 5) RETURN
648 IF(nstr(ig) == 0) THEN
649 lvers(ig)=0
650 ELSE
651 lvers(ig)=lvers(ig)+1
652 END IF
653 nstr(ig)=1 ! by GF
654 ! remove stored elements
655 IF(jflc(1,ig) /= 0) CALL stmarm(jflc(1,ig))
656 IF(kflc(1,ig) /= 0) CALL stmarm(kflc(1,ig))
657 igtp(ig)=ityp
658 gtext(ig)=text
659 DO j=1,5
660 jflc(j,ig)=0
661 END DO
662 jflc(5,ig)=nlimit
663 IF(ityp == 5) THEN
664 DO j=1,5
665 kflc(j,ig)=0
666 END DO
667 jflc(5,ig)=128 ! maximum of 128 values
668 kflc(5,ig)=narr
669 nst(1,ig)=0
670 nst(2,ig)=0
671 nst(3,ig)=1
672 DO j=1,10
673 xyplws(j,ig)=0.0
674 END DO
675 END IF
676 RETURN
677
678 entry gmpxy(ig,x,y) ! add (X,Y) pair
679 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG
680 IF(igtp(ig) < 1.OR.igtp(ig) > 3) RETURN ! check type
681 CALL stmapr(jflc(1,ig),x,y)
682 RETURN
683
684 entry gmpxyd(ig,x,y,dx,dy) ! add (X,Y,DX,DY)
685 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG
686 IF(igtp(ig) /= 4) RETURN
687 four(1)=x
688 four(2)=y
689 four(3)=dx
690 four(4)=dy
691 CALL stmadp(jflc(1,ig),four)
692 RETURN
693
694 entry gmpms(ig,x,y) ! mean sigma(X) from Y
695 ! mean sigma from Y, as a function of X
696 ! WRITE(*,*) 'GMPMS ',IG,X,Y
697
698 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG
699 IF(igtp(ig) /= 5) RETURN
700
701 xyplws(10,ig)=x ! last X coordinate
702 IF(nst(1,ig) == 0) THEN
703 y1=y
704 nst(1,ig)=1
705 IF(kflc(3,ig) == 0) xyplws(9,ig)=x ! start coordinate
706 ELSE
707 nst(1,ig)=0
708 CALL stmapr(kflc(1,ig),y1,y) ! store pair
709 IF(kflc(3,ig) >= kflc(5,ig)) THEN
710 CALL stmacp(kflc(1,ig),array,n) ! get data
711 CALL stmarm(kflc(1,ig)) ! remove data
712 n=n+n
713 CALL rmesig(array,n,xyplws(2,ig),xyplws(4,ig))
714 nst(2,ig)=nst(2,ig)+1
715 IF(nst(2,ig) == 1) xyplws(7,ig)=xyplws(9,ig)
716 xyplws(8,ig)=x ! end coordinate
717 xyplws(5,ig)=xyplws(5,ig)+xyplws(2,ig)
718 xyplws(6,ig)=xyplws(6,ig)+xyplws(4,ig)
719 IF(nst(2,ig) == nst(3,ig)) THEN
720 xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig))
721 xyplws(2,ig)=xyplws(5,ig)/real(nst(3,ig),mps)
722 xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig))
723 xyplws(4,ig)=xyplws(6,ig)/real(nst(3,ig),mps)
724 xyplws(5,ig)=0.0
725 xyplws(6,ig)=0.0
726 nst(2,ig)=0
727 CALL stmadp(jflc(1,ig),xyplws(1,ig))
728 IF(jflc(3,ig) >= jflc(5,ig)) THEN
729 CALL stmacp(jflc(1,ig),array4,n) ! get data
730 n=n/2
731 CALL stmarm(jflc(1,ig)) ! remove data
732 DO i=1,n,2 ! average
733 xyplws(7,ig)=array4(1,i )-array4(3, i)
734 xyplws(8,ig)=array4(1,i+1)+array4(3,i+1)
735 xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig))
736 xyplws(2,ig)=0.5*(array4(2,i)+array4(2,i+1))
737 xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig))
738 xyplws(4,ig)=0.5*(array4(4,i)+array4(4,i+1))
739 CALL stmadp(jflc(1,ig),xyplws(1,ig))
740 END DO
741 nst(3,ig)=nst(3,ig)+nst(3,ig)
742 END IF
743 END IF
744 END IF
745 END IF
746 RETURN
747
748 entry gmprnt(ig) ! print XY data
749 IF(ig == 0) THEN
750 iga=1
751 igb=numgxy
752 ELSE
753 IF(ig <= 0.OR.ig > numgxy) RETURN
754 iga=ig
755 igb=ig
756 END IF
757 DO igc=iga,igb
758
759 IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN
760 WRITE(*,*) ' '
761 WRITE(*,*) 'Store ',igc,': ',gtext(igc)
762 IF(jflc(4,igc) == 0) THEN
763 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
764 ELSE
765 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
766 jflc(3,igc),', ',jflc(4,igc)
767 END IF
768
769 CALL stmacp(jflc(1,igc),array,na) ! get all data
770
771 DO n=1,na
772 WRITE(*,102) n, array(1,n),array(2,n)
773 END DO
774
775 ELSE IF(igtp(igc) == 4) THEN
776
777 WRITE(*,*) ' '
778 WRITE(*,*) 'Store ',igc,': ',gtext(igc)
779 IF(jflc(4,igc) == 0) THEN
780 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
781 ELSE
782 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
783 jflc(3,igc),', ',jflc(4,igc)
784 END IF
785
786 CALL stmacp(jflc(1,igc),array,na) ! get all data
787 na=na/2
788
789 DO n=1,na
790 WRITE(*,102) n,(array4(j,n),j=1,4)
791 END DO
792
793 ELSE IF(igtp(igc) == 5) THEN
794
795 CALL stmacp(kflc(1,igc),array,n) ! get data
796 CALL stmarm(kflc(1,igc)) ! remove data
797 n=n+n
798 IF(nst(1,igc) == 1) THEN
799 n=n+1
800 array1(n)=y1
801 nst(1,igc)=0 ! reset
802 END IF
803 IF(n /= 0) THEN
804 xyplws(7,igc)=xyplws( 9,igc)
805 xyplws(8,igc)=xyplws(10,igc)
806 CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc))
807 wght=real(n,mps)/real(nst(3,igc)*kflc(5,igc),mps)
808 xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght
809 xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght
810 xyplws(2,igc)=xyplws(5,igc)/(real(nst(2,igc),mps)+wght)
811 xyplws(4,igc)=xyplws(6,igc)/(real(nst(2,igc),mps)+wght)
812 xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc))
813 xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc))
814 CALL stmadp(jflc(1,igc),xyplws(1,igc))
815 END IF
816
817 WRITE(*,*) ' '
818 WRITE(*,*) 'Store ',igc,': ',gtext(igc)
819 IF(jflc(4,igc) == 0) THEN
820 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc)
821 ELSE
822 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', &
823 jflc(3,igc),', ',jflc(4,igc)
824 END IF
825
826 CALL stmacp(jflc(1,igc),array,na) ! get all data
827 na=na/2
828 DO n=1,na
829 WRITE(*,102) n,(array4(j,n),j=1,4)
830 END DO
831 END IF
832 END DO
833 RETURN
834
835 entry gmplun(lunw) ! unit for output
836 lun=lunw
837 RETURN
838
839 entry gmpwrt(ig) ! write XY text file
840 IF(lun <= 0) RETURN
841 IF(ig == 0) THEN
842 iga=1
843 igb=numgxy
844 ELSE
845 IF(ig <= 0.OR.ig > numgxy) RETURN
846 iga=ig
847 igb=ig
848 END IF
849 DO igc=iga,igb
850 IF(igtp(igc) == 5) THEN
851
852 CALL stmacp(kflc(1,igc),array,n) ! get data
853 CALL stmarm(kflc(1,igc)) ! remove data
854 n=n+n
855 IF(nst(1,igc) == 1) THEN
856 n=n+1
857 array1(n)=y1
858 nst(1,igc)=0 ! reset
859 END IF
860 IF(n /= 0) THEN
861 xyplws(7,igc)=xyplws( 9,igc)
862 xyplws(8,igc)=xyplws(10,igc)
863 CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc))
864 wght=real(n,mps)/real(nst(3,igc)*kflc(5,igc),mps)
865 xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght
866 xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght
867 xyplws(2,igc)=xyplws(5,igc)/(real(nst(2,igc),mps)+wght)
868 xyplws(4,igc)=xyplws(6,igc)/(real(nst(2,igc),mps)+wght)
869 xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc))
870 xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc))
871 CALL stmadp(jflc(1,igc),xyplws(1,igc))
872 END IF
873
874 END IF
875 IF(jflc(3,igc)+jflc(4,igc) /= 0) THEN
876 WRITE(lun,204) ' '
877 WRITE(lun,201) igc,lvers(igc),igtp(igc)
878 WRITE(lun,204) gtext(igc)
879 WRITE(lun,203) jflc(3,igc)+jflc(4,igc)
880 CALL stmacp(jflc(1,igc),array,na) ! get all data
881 IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN
882 WRITE(lun,204) 'x-y'
883 DO n=1,na
884 WRITE(lun,205) array(1,n),array(2,n)
885 END DO
886 ELSE IF(igtp(igc) == 4.OR.igtp(igc) == 5) THEN
887 WRITE(lun,204) 'x-y-dx-dy'
888 na=na/2
889 DO n=1,na
890 WRITE(lun,205) (array4(j,n),j=1,4)
891 END DO
892 END IF
893 WRITE(lun,204) 'end of xy-data'
894 END IF
895 END DO
896102 FORMAT(i12,4g15.7)
897 ! 103 FORMAT(' Index ___X___ ___Y___ '/
898 ! + ' ----- -------------- --------------')
899 ! 104 FORMAT(' Index ___X___ ___Y___ ',
900 ! + ' ___DX__ ___DY__ '/
901 ! + ' ----- -------------- --------------',
902 ! + ' -------------- --------------')
903201 FORMAT('XY-Data ',i4,10x,'version ',i4,10x,'type',i2)
904203 FORMAT(10x,'stored not-stored ',2i10)
905204 FORMAT(a)
906205 FORMAT(3x,4g15.7)
907END SUBROUTINE gmpdef
908
909
910SUBROUTINE stmars ! init/reset storage
911 USE mpdef
912
913 IMPLICIT NONE
914 INTEGER(mpi) :: i
915 INTEGER(mpi) :: ifre
916 INTEGER(mpi) :: ifrea
917 INTEGER(mpi) :: ifreb
918 INTEGER(mpi) :: ind
919 INTEGER(mpi) :: j
920 INTEGER(mpi) :: n
921 REAL(mps) :: x
922 REAL(mps) :: y
923 INTEGER(mpi), PARAMETER :: ndim=5000 ! storage dimension, should be NUMGXY*NLIMIT
924 REAL(mps) :: tk(2,ndim) ! pair storage for data pairs
925 INTEGER(mpi) :: next(ndim) ! pointer
926 INTEGER(mpi) :: iflc1 ! first and last index of free pairs
927 INTEGER(mpi) ::iflc2 ! first and last index of free pairs
928 SAVE
929
930 REAL(mps) :: four(4) ! double_pair, copy array
931 REAL(mps) ::array(2,*) ! double_pair, copy array
932 INTEGER(mpi) :: jflc(5) ! user array
933 ! JFLC(1) = first used index
934 ! JFLC(2) = last used index
935 ! JFLC(3) = counter of used places
936 ! JFLC(4) = counter of ignored
937 ! JFLC(5) = limit for JFLC(3)
938 ! ...
939 DO i=1,ndim
940 next(i)=i+1 ! pointer to next free location
941 tk(1,i)=0.0 ! reset
942 tk(2,i)=0.0
943 END DO
944 next(ndim)=0 ! ... and end pointer
945 iflc1=1 ! index first free pair
946 iflc2=ndim ! index last free pair
947 RETURN
948
949 entry stmapr(jflc,x,y) ! store pair (X,Y)
950 ifre=iflc1 ! index of free place
951 IF(ifre == 0.OR.jflc(3) >= jflc(5)) THEN ! overflow
952 jflc(4)=jflc(4)+1
953 ELSE
954 iflc1=next(ifre) ! pointer to new free location
955 IF(jflc(1) == 0) THEN ! first item
956 jflc(1)=ifre
957 ELSE
958 next(jflc(2))=ifre
959 END IF
960 next(ifre)=0
961 jflc(2)=ifre ! last index
962 jflc(3)=jflc(3)+1 ! counter
963 tk(1,ifre)=x
964 tk(2,ifre)=y
965 END IF
966 RETURN
967
968 entry stmadp(jflc,four) ! store double pair
969 ifrea=iflc1 ! index of 1. free place
970 IF(ifrea == 0) THEN ! overflow
971 jflc(4)=jflc(4)+1
972 ELSE
973 ifreb=next(iflc1) ! index of 2. free place
974 IF(ifreb == 0.OR.jflc(3) >= 2*jflc(5)) THEN ! overflow
975 jflc(4)=jflc(4)+1
976 ELSE
977 iflc1=next(ifreb) ! pointer to new free location
978 IF(jflc(1) == 0) THEN ! first item
979 jflc(1)=ifrea
980 ELSE
981 next(jflc(2))=ifrea
982 END IF
983 next(ifreb)=0
984 jflc(2)=ifreb ! last index
985 jflc(3)=jflc(3)+1 ! counter
986 tk(1,ifrea)=four(1)
987 tk(2,ifrea)=four(2)
988 tk(1,ifreb)=four(3)
989 tk(2,ifreb)=four(4)
990 END IF
991 END IF
992 RETURN
993
994 entry stmacp(jflc,array,n) ! copy (cp) all pairs to array
995 n=0
996 ind=jflc(1)
99710 IF(ind == 0) RETURN
998 n=n+1
999 array(1,n)=tk(1,ind)
1000 array(2,n)=tk(2,ind)
1001 ind=next(ind)
1002 GO TO 10
1003
1004 entry stmarm(jflc) ! remove (rm) stored paiirs
1005 next(iflc2)=jflc(1) ! connect to free space
1006 iflc2=jflc(2) ! new last free index
1007 DO j=1,4
1008 jflc(j)=0
1009 END DO
1010END SUBROUTINE stmars ! init/
1011
1012SUBROUTINE rmesig(x,n,xloc,xsca) ! robust mean and sigma
1013 USE mpdef
1014
1015 IMPLICIT NONE
1016 INTEGER(mpi) :: i
1017 ! robust determination of location and scale parameter,
1018 ! for Gaussian data: location=mean and scale=standard deviation
1019 ! XLOC = median of X_i (N values in array X(N))
1020 ! XCSA = median of | X_i - XLOC |, times 1.4826
1021
1022 REAL(mps), INTENT(IN OUT) :: x(n) ! input array, modified
1023 INTEGER(mpi), INTENT(IN) :: n
1024 REAL(mps), INTENT(OUT) :: xloc
1025 REAL(mps), INTENT(OUT) :: xsca
1026 SAVE
1027 ! ...
1028 xloc=0.0
1029 xsca=0.0
1030 IF(n <= 0) RETURN
1031 CALL heapf(x,n) ! sort
1032 xloc=0.5*(x((n+1)/2)+x((n+2)/2)) ! location
1033 DO i=1,n
1034 x(i)=abs(x(i)-xloc)
1035 END DO
1036 CALL heapf(x,n) ! sort
1037 xsca=1.4826*0.5*(x((n+1)/2)+x((n+2)/2)) ! dispersion
1038END SUBROUTINE rmesig
1039
1040
1041
hmpmak
subroutine hmpmak(inhist, fnhist, jnhist, xl, dl)
Definition mphistab.f90:351
stmars
subroutine stmars
Definition mphistab.f90:911
gmpdef
subroutine gmpdef(ig, ityp, text)
Definition mphistab.f90:582
bintab
subroutine bintab(tab, n, xa, xb)
Definition mphistab.f90:416
hmpdef
subroutine hmpdef(ih, xa, xb, text)
Definition mphistab.f90:93
rmesig
subroutine rmesig(x, n, xloc, xsca)
Definition mphistab.f90:1013
kprint
subroutine kprint(lun, list, n)
Definition mphistab.f90:520
heapf
subroutine heapf(a, n)
Heap sort direct (real).
Definition mpnum.f90:1468
mpdef
Definition of constants.
Definition mpdef.f90:24
mpdef::mpi
integer, parameter mpi
Definition mpdef.f90:34
pfvert
subroutine pfvert(n, x)
Vertical print of floating point data.
Definition vertpr.f90:225
pivert
subroutine pivert(n, list)
Vertical print of integer data.
Definition vertpr.f90:185
psvert
subroutine psvert(xa, xb)
Print scale.
Definition vertpr.f90:259