 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:07:45 00004 00070 00071 ! *************************** Histograms ****************************** 00072 00073 SUBROUTINE hmpdef(ih,xa,xb,text) ! book, reset histogram 00074 IMPLICIT NONE 00075 INTEGER:: i 00076 INTEGER:: iha 00077 INTEGER:: ihb 00078 INTEGER:: ihc 00079 INTEGER:: ix 00080 INTEGER:: j 00081 INTEGER:: lun 00082 INTEGER:: lunw 00083 INTEGER:: nbin 00084 INTEGER:: nn 00085 REAL:: x 00086 REAL:: xcent 00087 REAL:: xmean 00088 REAL:: xsigm 00089 ! book millepede histogram, 120 bins 00090 00091 INTEGER, INTENT(IN) :: ih 00092 REAL, INTENT(IN) :: xa 00093 REAL, INTENT(IN) :: xb 00094 CHARACTER (LEN=*), INTENT(IN) :: text 00095 INTEGER, PARAMETER :: numhis=15 00096 INTEGER:: inhist(120,numhis) 00097 INTEGER::jnhist(5,numhis) 00098 INTEGER::khist(numhis) 00099 REAL:: fnhist(120,numhis) 00100 EQUIVALENCE (inhist(1,1),fnhist(1,1)) 00101 INTEGER:: kvers(numhis) 00102 REAL:: xl(6,numhis) 00103 DOUBLE PRECISION:: dl(2,numhis) 00104 CHARACTER (LEN=60):: htext(numhis) 00105 SAVE 00106 DATA khist/numhis*0/,lun/7/ 00107 ! ... 00108 IF(ih <= 0.OR.ih > numhis) RETURN 00109 ! IF(XA.EQ.XB) RETURN 00110 DO i=1,120 00111 inhist(i,ih)=0 00112 END DO 00113 DO j=1,5 00114 jnhist(j,ih)=0 00115 END DO 00116 xl(1,ih)=xa 00117 xl(2,ih)=xb 00118 xl(3,ih)=0.0 00119 IF(xa /= xb) xl(3,ih)=120.0/(xb-xa) 00120 xl(6,ih)=0.5*(xa+xb) ! center 00121 IF(khist(ih) == 0) THEN 00122 kvers(ih)=0 00123 ELSE 00124 kvers(ih)=kvers(ih)+1 00125 END IF 00126 khist(ih)=1 ! flt.pt. (lin) 00127 htext(ih)=text 00128 dl(1,ih)=0.0D0 00129 dl(2,ih)=0.0D0 00130 RETURN 00131 00132 ENTRY hmpldf(ih,text) ! book, reset log histogram 00133 IF(ih <= 0.OR.ih > numhis) RETURN 00134 DO i=1,120 00135 inhist(i,ih)=0 00136 END DO 00137 DO j=1,5 00138 jnhist(j,ih)=0 00139 END DO 00140 IF(khist(ih) == 0) THEN 00141 kvers(ih)=0 00142 ELSE 00143 kvers(ih)=kvers(ih)+1 00144 END IF 00145 khist(ih)=2 ! integer log 00146 htext(ih)=text 00147 xl(1,ih)=0.0 00148 xl(2,ih)=6.0 00149 RETURN 00150 00151 ENTRY hmpent(ih,x) ! entry flt.pt. 00152 IF(ih <= 0.OR.ih > numhis) RETURN 00153 IF(khist(ih) /= 1) RETURN 00154 IF(jnhist(4,ih) >= 2147483647) RETURN 00155 jnhist(4,ih)=jnhist(4,ih)+1 ! count 00156 IF(jnhist(4,ih) <= 120) THEN 00157 fnhist(jnhist(4,ih),ih)=x ! store value 00158 IF(jnhist(4,ih) == 120) THEN 00159 CALL hmpmak(inhist(1,ih),fnhist(1,ih),jnhist(1,ih), xl(1,ih),dl(1,ih)) 00160 END IF 00161 RETURN 00162 END IF 00163 ! IF(JNHIST(1,IH)+JNHIST(2,IH)+JNHIST(3,IH).EQ.0) THEN 00164 ! XL(4,IH)=X 00165 ! XL(5,IH)=X 00166 ! END IF 00167 i=INT(1.0+xl(3,ih)*(x-xl(1,ih))) ! X - Xmin 00168 j=2 00169 IF(i < 1) j=1 00170 IF(i > 120) j=3 00171 jnhist(j,ih)=jnhist(j,ih)+1 00172 xl(4,ih)=MIN(xl(4,ih),x) 00173 xl(5,ih)=MAX(xl(5,ih),x) 00174 IF(j /= 2) RETURN 00175 inhist(i,ih)=inhist(i,ih)+1 00176 dl(1,ih)=dl(1,ih)+ x-xl(6,ih) 00177 dl(2,ih)=dl(2,ih)+(x-xl(6,ih))**2 00178 RETURN 00179 00180 ENTRY hmplnt(ih,ix) ! entry integer 00181 IF(ih <= 0.OR.ih > numhis) RETURN 00182 IF(khist(ih) /= 2) RETURN 00183 IF(jnhist(1,ih) >= 2147483647) RETURN 00184 IF(ix <= 0) THEN 00185 jnhist(1,ih)=jnhist(1,ih)+1 00186 ELSE 00187 IF(jnhist(4,ih) == 0) jnhist(4,ih)=ix 00188 IF(jnhist(5,ih) == 0) jnhist(5,ih)=ix 00189 jnhist(4,ih)=MIN(jnhist(4,ih),ix) 00190 jnhist(5,ih)=MAX(jnhist(5,ih),ix) 00191 i=INT(1.0+20.0*LOG10(REAL(ix))) 00192 j=2 00193 IF(i < 1) j=1 00194 IF(i > 120) j=3 00195 IF(j == 2) inhist(i,ih)=inhist(i,ih)+1 00196 jnhist(j,ih)=jnhist(j,ih)+1 00197 END IF 00198 RETURN 00199 00200 ENTRY hmprnt(ih) ! print, content vert 00201 IF(ih == 0) THEN 00202 iha=1 00203 ihb=numhis 00204 ELSE 00205 IF(ih <= 0.OR.ih > numhis) RETURN 00206 iha=ih 00207 ihb=ih 00208 END IF 00209 DO ihc=iha,ihb 00210 IF(khist(ihc) /= 0) THEN 00211 IF(khist(ihc) == 1) THEN 00212 CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), & 00213 xl(1,ihc),dl(1,ihc)) 00214 END IF 00215 nn=jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc) 00216 IF(nn /= 0.OR.khist(ihc) == 3) THEN 00217 WRITE(*,111) 00218 111 FORMAT(' ______',2('______________________________')) 00219 IF(kvers(ihc) == 1) THEN 00220 WRITE(*,*) 'Histogram',ihc,': ',htext(ihc) 00221 ELSE 00222 WRITE(*,*) 'Histogram',ihc,'/',kvers(ihc),': ',htext(ihc) 00223 END IF 00224 IF(khist(ihc) == 1) THEN 00225 WRITE(*,*) ' Out_low inside out_high = ', (jnhist(j,ihc),j=1,3) 00226 ELSE IF(khist(ihc) == 2) THEN 00227 WRITE(*,*) ' 0_or_negative inside above_10^6 = ', & 00228 (jnhist(j,ihc),j=1,3) 00229 END IF 00230 IF(khist(ihc) == 3) THEN 00231 CALL pfvert(120,fnhist(1,ihc)) 00232 END IF 00233 IF(jnhist(2,ihc) /= 0) THEN ! integer content 00234 CALL pivert(120,inhist(1,ihc)) 00235 IF(khist(ihc) == 1) THEN 00236 CALL psvert(xl(1,ihc),xl(2,ihc)) 00237 ELSE IF(khist(ihc) == 2) THEN 00238 CALL psvert(0.0,6.0) 00239 END IF 00240 END IF 00241 IF(khist(ihc) == 1) THEN 00242 WRITE(*,*) ' Min and Max are',xl(4,ihc),xl(5,ihc) 00243 IF(jnhist(2,ihc) > 1) THEN 00244 xmean=REAL(xl(6,ihc)+dl(1,ihc)/FLOAT(jnhist(2,ihc))) 00245 xcent=0.5*(xl(1,ihc)+xl(2,ihc)) 00246 xsigm=REAL((dl(2,ihc)-dl(1,ihc)**2/FLOAT(jnhist(2,ihc)))) 00247 xsigm=SQRT(xsigm/FLOAT(jnhist(2,ihc)-1)) 00248 WRITE(*,*) ' Mean and sigma are', xmean,' +-',xsigm 00249 END IF 00250 ELSE IF(khist(ihc) == 2) THEN 00251 WRITE(*,*) ' Plot of log10 of entries. Min and Max are', & 00252 jnhist(4,ihc),jnhist(5,ihc) 00253 END IF 00254 END IF 00255 END IF 00256 END DO 00257 RETURN 00258 00259 ENTRY hmplun(lunw) ! unit for output 00260 lun=lunw 00261 RETURN 00262 00263 ENTRY hmpwrt(ih) ! write histogram text file 00264 IF(lun <= 0) RETURN 00265 IF(ih == 0) THEN 00266 iha=1 00267 ihb=numhis 00268 ELSE 00269 IF(ih <= 0.OR.ih > numhis) RETURN 00270 iha=ih 00271 ihb=ih 00272 END IF 00273 00274 DO ihc=iha,ihb ! histogram loop 00275 IF(khist(ihc) /= 0) THEN 00276 IF(khist(ihc) == 1) THEN 00277 CALL hmpmak(inhist(1,ihc),fnhist(1,ihc),jnhist(1,ihc), & 00278 xl(1,ihc),dl(1,ihc)) 00279 END IF 00280 nbin=120 00281 WRITE(lun,204) ' ' 00282 WRITE(lun,201) ihc,kvers(ihc),khist(ihc) 00283 WRITE(lun,204) htext(ihc) 00284 IF (jnhist(1,ihc)+jnhist(2,ihc)+jnhist(3,ihc) == 0 & 00285 .AND.xl(1,ihc) == xl(2,ihc)) THEN 00286 ! hist is empty and hist range makes no sense 00287 ! - cause: hist with 'variable edges' was never filled 00288 ! - workaround: make lower and upper edge of hist differ in output 00289 WRITE(lun,202) nbin,xl(1,ihc)-0.001,xl(2,ihc)+0.001 00290 ELSE 00291 WRITE(lun,202) nbin,xl(1,ihc),xl(2,ihc) 00292 END IF 00293 WRITE(lun,203) (jnhist(j,ihc),j=1,3) 00294 WRITE(lun,204) 'bincontent' 00295 IF(khist(ihc) == 1.OR.khist(ihc) == 2) THEN 00296 CALL kprint(lun,inhist(1,ihc),nbin) 00297 ELSE 00298 WRITE(lun,219) (fnhist(i,ihc),i=1,nbin) 00299 END IF 00300 00301 IF(khist(ihc) == 1) THEN 00302 WRITE(lun,205) xl(4,ihc),xl(5,ihc) 00303 ELSE IF(khist(ihc) == 2) THEN 00304 WRITE(lun,205) FLOAT(jnhist(4,ihc)),FLOAT(jnhist(5,ihc)) 00305 END IF 00306 IF(khist(ihc) == 1) THEN 00307 IF(jnhist(2,ihc) > 1) THEN 00308 xmean=REAL(xl(6,ihc)+dl(1,ihc)/FLOAT(jnhist(2,ihc))) 00309 xcent=0.5*(xl(1,ihc)+xl(2,ihc)) 00310 xsigm=REAL((dl(2,ihc)-dl(1,ihc)**2/FLOAT(jnhist(2,ihc)))) 00311 xsigm=SQRT(xsigm/FLOAT(jnhist(2,ihc)-1)) 00312 WRITE(lun,206) xmean,xsigm 00313 END IF 00314 END IF 00315 WRITE(lun,204) 'end of histogram' 00316 END IF 00317 END DO 00318 00319 201 FORMAT('Histogram ',i4,10X,'version ',i4,10X,'type',i2) 00320 202 FORMAT(10X,' bins, limits ',i4,2G15.5) 00321 203 FORMAT(10X,'out-low inside out-high ',3I10) 00322 204 FORMAT(a) 00323 205 FORMAT('minmax',2E15.7) 00324 206 FORMAT('meansigma',2E15.7) 00325 00326 219 FORMAT(4E15.7) 00327 END SUBROUTINE hmpdef 00328 00329 SUBROUTINE hmpmak(inhist,fnhist,jnhist,xl,dl) ! hist scale from data 00330 IMPLICIT NONE 00331 INTEGER:: i 00332 INTEGER:: j 00333 INTEGER:: k 00334 INTEGER:: nn 00335 REAL:: x 00336 REAL:: xa 00337 REAL:: xb 00338 00339 INTEGER, INTENT(OUT) :: inhist(120) 00340 REAL, INTENT(IN) :: fnhist(120) 00341 INTEGER, INTENT(IN OUT) :: jnhist(5) 00342 REAL, INTENT(IN OUT) :: xl(6) 00343 DOUBLE PRECISION, INTENT(OUT) :: dl(2) 00344 REAL:: cphist(120) 00345 00346 00347 00348 SAVE 00349 ! ... 00350 nn=jnhist(4) 00351 ! WRITE(*,*) 'HMPMAK: NN,JNHIST(5)',NN,JNHIST(5) 00352 IF(nn == 0.OR.jnhist(5) /= 0) RETURN 00353 jnhist(5)=1 00354 DO i=1,nn 00355 ! WRITE(*,*) 'copy ',I,FNHIST(I) 00356 cphist(i)=fnhist(i) 00357 END DO 00358 CALL heapf(cphist,nn) 00359 IF(xl(3) == 0.0) THEN 00360 CALL bintab(cphist,nn,xa,xb) 00361 xl(1)=xa 00362 xl(2)=xb 00363 xl(3)=0.0 00364 IF(xa /= xb) xl(3)=120.0/(xb-xa) 00365 xl(6)=0.5*(xa+xb) ! center 00366 END IF 00367 xl(4)=cphist( 1) 00368 xl(5)=cphist(nn) 00369 ! WRITE(*,*) 'XL ',XL 00370 DO i=1,nn 00371 inhist(i)=0 00372 END DO 00373 DO k=1,nn 00374 x=cphist(k) 00375 i=INT(1.0+xl(3)*(x-xl(1))) ! X - Xmin 00376 ! WRITE(*,*) 'K,I,X ',K,I,X 00377 j=2 00378 IF(i < 1) j=1 00379 IF(i > 120) j=3 00380 jnhist(j)=jnhist(j)+1 00381 IF(j == 2) THEN 00382 inhist(i)=inhist(i)+1 00383 dl(1)=dl(1)+ x-xl(6) 00384 dl(2)=dl(2)+(x-xl(6))**2 00385 END IF 00386 END DO 00387 END SUBROUTINE hmpmak 00388 00389 00390 00391 00392 SUBROUTINE bintab(tab,n,xa,xb) ! hist scale from data 00393 IMPLICIT NONE 00394 REAL:: dd 00395 REAL:: dx 00396 INTEGER:: i 00397 INTEGER:: iexp 00398 INTEGER:: ii 00399 INTEGER:: j 00400 INTEGER:: m1 00401 INTEGER:: m2 00402 INTEGER:: n1 00403 INTEGER:: n2 00404 REAL:: rat 00405 REAL:: x1 00406 REAL:: x2 00407 REAL:: xx 00408 ! Bin limits XA and XB from TAB(N) 00409 00410 REAL, INTENT(IN) :: tab(n) 00411 INTEGER, INTENT(IN) :: n 00412 REAL, INTENT(OUT) :: xa 00413 REAL, INTENT(OUT) :: xb 00414 00415 REAL:: bin(10) 00416 DATA bin/1.0,1.5,2.0,3.0,4.0,5.0,8.0,10.0,15.0,20.0/ 00417 SAVE 00418 ! ... 00419 00420 CALL heapf(tab,n) ! reduced statistic 00421 ! WRITE(*,*) ' ' 00422 ! WRITE(*,*) 'Sorted ',(TAB(I),I=1,N) 00423 IF(n < 100) THEN 00424 x1=tab(1) 00425 x2=tab(n) 00426 ! WRITE(*,*) 'reduced statistic X1 X2 ',X1,X2 00427 ELSE ! large statistic 00428 m1=INT(1.0+0.05*FLOAT(n)) 00429 m2=INT(1.0+0.16*FLOAT(n)) 00430 x1=tab(m1)-4.0*(tab(m2)-tab(m1)) 00431 IF(x1 < 0.0.AND.tab(1) >= 0.0) x1=tab(1) 00432 x2=tab(n+1-m1)+4.0*(tab(n+1-m1)-tab(n+1-m2)) 00433 IF(x2 > 0.0.AND.tab(n) <= 0.0) x2=tab(n) 00434 ! WRITE(*,*) 'large statistic ',X1,X2 00435 ! WRITE(*,*) 'min und max ',TAB(1),TAB(N) 00436 IF(x1*tab(1) <= 0.0) x1=0.0 00437 IF(x2*tab(n) <= 0.0) x2=0.0 00438 ! WRITE(*,*) 'large statistic zero ',X1,X2 00439 IF(x1*x2 < 0.0.AND.MIN(-x1,x2) > 0.6*MAX(-x1,x2)) THEN 00440 xx=MAX(-x1,x2) ! symmetry 00441 x1=-xx 00442 x2=+xx 00443 ELSE IF(x1*x2 > 0.0.AND. & ! include zero ? 00444 ABS(MIN(x1,x2)) < 0.4*ABS(MAX(x1,x2))) THEN 00445 IF(x1 < 0.0) THEN 00446 x2=0.0 00447 ELSE 00448 x1=0.0 00449 END IF 00450 END IF 00451 ! WRITE(*,*) 'large statistic ',X1,X2 00452 END IF 00453 IF(x1 == x2) THEN 00454 x1=x1-1.0 00455 x2=x2+1.0 00456 END IF 00457 dx=x2-x1 00458 ! WRITE(*,*) 'X1,X2,DX ',X1,X2,DX 00459 rat=0.0 00460 ii=1 00461 DO j=1,11 00462 i=j 00463 IF(j == 11) i=ii 00464 iexp=INT(101.0+LOG10(dx)-LOG10(6.0*bin(i))) 00465 iexp=iexp-100 00466 dd=bin(i)*10.0**iexp 00467 00468 n1=INT(ABS(x1)/dd) 00469 IF(x1 < 0.0) n1=-n1 00470 IF(FLOAT(n1)*dd > x1) n1=n1-1 00471 ! WRITE(*,*) 'Bin ',I,N1,N1*DD,X1 00472 00473 n2=INT(ABS(x2)/dd) 00474 IF(x2 < 0.0) n2=-n2 00475 IF(FLOAT(n2)*dd < x2) n2=n2+1 00476 ! WRITE(*,*) 'Bin ',I,N2,N2*DD,X2 00477 10 IF(n2-n1 < 6) THEN 00478 IF(n1 /= 0) n1=n1-1 00479 IF(n2-n1 < 6.AND.n2 /= 0) n2=n2+1 00480 GO TO 10 00481 END IF 00482 ! WRITE(*,*) 'corrected N1 N2 ',N1,N2 00483 xa=SIGN(FLOAT(n1)*dd,x1) 00484 xb=SIGN(FLOAT(n2)*dd,x2) 00485 ! WRITE(*,*) J,' resulting limits XA XB ',XA,XB 00486 IF((x2-x1)/(xb-xa) > rat) THEN 00487 ii=i 00488 rat=(x2-x1)/(xb-xa) 00489 END IF 00490 END DO 00491 ! WRITE(*,*) J,' resulting limits XA XB ',XA,XB 00492 END SUBROUTINE bintab 00493 00494 SUBROUTINE kprint(lun,list,n) ! print integer array 00495 IMPLICIT NONE 00496 INTEGER:: i 00497 INTEGER:: ia 00498 INTEGER:: ib 00499 INTEGER:: k 00500 INTEGER:: ln 00501 INTEGER:: lp 00502 INTEGER:: np 00503 ! print integer array LIST(N) 00504 00505 INTEGER, INTENT(IN OUT) :: lun 00506 INTEGER, INTENT(IN) :: list(n) 00507 INTEGER, INTENT(IN) :: n 00508 INTEGER:: li(7) 00509 DATA li/2,3,4,6,8,9,12/ ! number of characters 00510 SAVE 00511 ! ... 00512 ib=0 00513 10 ia=ib+1 00514 IF(ia > n) RETURN 00515 DO k=1,7 00516 np=72/li(k) 00517 ib=MIN(ia-1+np,n) 00518 IF(k <= 6) THEN 00519 lp=10**(li(k)-1)-1 ! maximum positive 00520 ln=-lp/10 ! minimum negative 00521 DO i=ia,ib 00522 IF(list(i) > lp.OR.list(i) < ln) GO TO 20 00523 END DO 00524 END IF 00525 IF(k == 1) THEN 00526 WRITE(lun,101) (list(i),i=ia,ib) 00527 ELSE IF(k == 2) THEN 00528 WRITE(lun,102) (list(i),i=ia,ib) 00529 ELSE IF(k == 3) THEN 00530 WRITE(lun,103) (list(i),i=ia,ib) 00531 ELSE IF(k == 4) THEN 00532 WRITE(lun,104) (list(i),i=ia,ib) 00533 ELSE IF(k == 5) THEN 00534 WRITE(lun,105) (list(i),i=ia,ib) 00535 ELSE IF(k == 6) THEN 00536 WRITE(lun,106) (list(i),i=ia,ib) 00537 ELSE IF(k == 7) THEN 00538 WRITE(lun,107) (list(i),i=ia,ib) 00539 END IF 00540 GO TO 10 00541 20 CONTINUE 00542 END DO 00543 101 FORMAT(36I2) 00544 102 FORMAT(24I3) 00545 103 FORMAT(18I4) 00546 104 FORMAT(12I6) 00547 105 FORMAT( 9I8) 00548 106 FORMAT( 8I9) 00549 107 FORMAT( 6I12) 00550 END SUBROUTINE kprint 00551 00552 ! ***************************** XY data **************************** 00553 00554 SUBROUTINE gmpdef(ig,ityp,text) ! book, reset XY storage 00555 IMPLICIT NONE 00556 REAL:: dx 00557 REAL:: dy 00558 INTEGER:: i 00559 INTEGER:: iga 00560 INTEGER:: igb 00561 INTEGER:: igc 00562 INTEGER:: j 00563 INTEGER:: lun 00564 INTEGER:: lunw 00565 INTEGER:: n 00566 INTEGER:: na 00567 REAL:: wght 00568 REAL:: x 00569 REAL:: y 00570 REAL:: y1 00571 ! ITYP = 1 X,Y as dots 00572 ! = 2 X,Y as line 00573 ! = 3 X,Y as line and dots 00574 ! = 4 X,Y, DX,DY symbols 00575 00576 INTEGER, INTENT(IN) :: ig 00577 INTEGER, INTENT(IN) :: ityp 00578 CHARACTER (LEN=*), INTENT(IN) :: text 00579 INTEGER, PARAMETER :: narr=1000 00580 REAL:: array(2,narr) 00581 REAL::array4(4,narr/2) 00582 REAL::array1(narr+narr) 00583 REAL::four(4) 00584 EQUIVALENCE (array(1,1),array4(1,1),array1(1)) 00585 INTEGER, PARAMETER :: numgxy=10 00586 INTEGER, PARAMETER :: nlimit=500 00587 INTEGER:: nstr(numgxy) 00588 INTEGER::igtp(numgxy) 00589 INTEGER::lvers(numgxy) 00590 INTEGER::nst(3,numgxy) 00591 REAL:: xyplws(10,numgxy) 00592 INTEGER:: jflc(5,numgxy) 00593 INTEGER::kflc(5,numgxy) 00594 ! JFLC(1,.) = first used index 00595 ! JFLC(2,.) = last used index 00596 ! JFLC(3,.) = counter of used places 00597 ! JFLC(4,.) = counter of ignored 00598 ! JFLC(5,.) = limit for JFLC(3) 00599 CHARACTER (LEN=60):: gtext(numgxy) 00600 00601 LOGICAL:: start 00602 SAVE 00603 DATA start/.TRUE./,lun/7/ 00604 DATA nstr/numgxy*0/ ! by GF 00605 ! ... 00606 IF(start) THEN 00607 start=.FALSE. 00608 CALL stmars ! initialize storage 00609 DO i=1,numgxy 00610 DO j=1,5 00611 jflc(j,i)=0 00612 kflc(j,i)=0 00613 END DO 00614 END DO 00615 END IF 00616 00617 IF(ig < 1.OR.ig > numgxy) RETURN 00618 IF(ityp < 1.OR.ityp > 5) RETURN 00619 IF(nstr(ig) == 0) THEN 00620 lvers(ig)=0 00621 ELSE 00622 lvers(ig)=lvers(ig)+1 00623 END IF 00624 nstr(ig)=1 ! by GF 00625 ! remove stored elements 00626 IF(jflc(1,ig) /= 0) CALL stmarm(jflc(1,ig)) 00627 IF(kflc(1,ig) /= 0) CALL stmarm(kflc(1,ig)) 00628 igtp(ig)=ityp 00629 gtext(ig)=text 00630 DO j=1,5 00631 jflc(j,ig)=0 00632 END DO 00633 jflc(5,ig)=nlimit 00634 IF(ityp == 5) THEN 00635 DO j=1,5 00636 kflc(j,ig)=0 00637 END DO 00638 jflc(5,ig)=128 ! maximum of 128 values 00639 kflc(5,ig)=narr 00640 nst(1,ig)=0 00641 nst(2,ig)=0 00642 nst(3,ig)=1 00643 DO j=1,10 00644 xyplws(j,ig)=0.0 00645 END DO 00646 END IF 00647 RETURN 00648 00649 ENTRY gmpxy(ig,x,y) ! add (X,Y) pair 00650 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG 00651 IF(igtp(ig) < 1.OR.igtp(ig) > 3) RETURN ! check type 00652 CALL stmapr(jflc(1,ig),x,y) 00653 RETURN 00654 00655 ENTRY gmpxyd(ig,x,y,dx,dy) ! add (X,Y,DX,DY) 00656 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG 00657 IF(igtp(ig) /= 4) RETURN 00658 four(1)=x 00659 four(2)=y 00660 four(3)=dx 00661 four(4)=dy 00662 CALL stmadp(jflc(1,ig),four) 00663 RETURN 00664 00665 ENTRY gmpms(ig,x,y) ! mean sigma(X) from Y 00666 ! mean sigma from Y, as a function of X 00667 ! WRITE(*,*) 'GMPMS ',IG,X,Y 00668 00669 IF(ig < 1.OR.ig > numgxy) RETURN ! check argument IG 00670 IF(igtp(ig) /= 5) RETURN 00671 00672 xyplws(10,ig)=x ! last X coordinate 00673 IF(nst(1,ig) == 0) THEN 00674 y1=y 00675 nst(1,ig)=1 00676 IF(kflc(3,ig) == 0) xyplws(9,ig)=x ! start coordinate 00677 ELSE 00678 nst(1,ig)=0 00679 CALL stmapr(kflc(1,ig),y1,y) ! store pair 00680 IF(kflc(3,ig) >= kflc(5,ig)) THEN 00681 CALL stmacp(kflc(1,ig),array,n) ! get data 00682 CALL stmarm(kflc(1,ig)) ! remove data 00683 n=n+n 00684 CALL rmesig(array,n,xyplws(2,ig),xyplws(4,ig)) 00685 nst(2,ig)=nst(2,ig)+1 00686 IF(nst(2,ig) == 1) xyplws(7,ig)=xyplws(9,ig) 00687 xyplws(8,ig)=x ! end coordinate 00688 xyplws(5,ig)=xyplws(5,ig)+xyplws(2,ig) 00689 xyplws(6,ig)=xyplws(6,ig)+xyplws(4,ig) 00690 IF(nst(2,ig) == nst(3,ig)) THEN 00691 xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig)) 00692 xyplws(2,ig)=xyplws(5,ig)/FLOAT(nst(3,ig)) 00693 xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig)) 00694 xyplws(4,ig)=xyplws(6,ig)/FLOAT(nst(3,ig)) 00695 xyplws(5,ig)=0.0 00696 xyplws(6,ig)=0.0 00697 nst(2,ig)=0 00698 CALL stmadp(jflc(1,ig),xyplws(1,ig)) 00699 IF(jflc(3,ig) >= jflc(5,ig)) THEN 00700 CALL stmacp(jflc(1,ig),array4,n) ! get data 00701 n=n/2 00702 CALL stmarm(jflc(1,ig)) ! remove data 00703 DO i=1,n,2 ! average 00704 xyplws(7,ig)=array4(1,i )-array4(3, i) 00705 xyplws(8,ig)=array4(1,i+1)+array4(3,i+1) 00706 xyplws(1,ig)=0.5*(xyplws(7,ig)+xyplws(8,ig)) 00707 xyplws(2,ig)=0.5*(array4(2,i)+array4(2,i+1)) 00708 xyplws(3,ig)=0.5*(xyplws(8,ig)-xyplws(7,ig)) 00709 xyplws(4,ig)=0.5*(array4(4,i)+array4(4,i+1)) 00710 CALL stmadp(jflc(1,ig),xyplws(1,ig)) 00711 END DO 00712 nst(3,ig)=nst(3,ig)+nst(3,ig) 00713 END IF 00714 END IF 00715 END IF 00716 END IF 00717 RETURN 00718 00719 ENTRY gmprnt(ig) ! print XY data 00720 IF(ig == 0) THEN 00721 iga=1 00722 igb=numgxy 00723 ELSE 00724 IF(ig <= 0.OR.ig > numgxy) RETURN 00725 iga=ig 00726 igb=ig 00727 END IF 00728 DO igc=iga,igb 00729 00730 IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN 00731 WRITE(*,*) ' ' 00732 WRITE(*,*) 'Store ',igc,': ',gtext(igc) 00733 IF(jflc(4,igc) == 0) THEN 00734 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc) 00735 ELSE 00736 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', & 00737 jflc(3,igc),', ',jflc(4,igc) 00738 END IF 00739 00740 CALL stmacp(jflc(1,igc),array,na) ! get all data 00741 00742 DO n=1,na 00743 WRITE(*,102) n, array(1,n),array(2,n) 00744 END DO 00745 00746 ELSE IF(igtp(igc) == 4) THEN 00747 00748 WRITE(*,*) ' ' 00749 WRITE(*,*) 'Store ',igc,': ',gtext(igc) 00750 IF(jflc(4,igc) == 0) THEN 00751 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc) 00752 ELSE 00753 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', & 00754 jflc(3,igc),', ',jflc(4,igc) 00755 END IF 00756 00757 CALL stmacp(jflc(1,igc),array,na) ! get all data 00758 na=na/2 00759 00760 DO n=1,na 00761 WRITE(*,102) n,(array4(j,n),j=1,4) 00762 END DO 00763 00764 ELSE IF(igtp(igc) == 5) THEN 00765 00766 CALL stmacp(kflc(1,igc),array,n) ! get data 00767 CALL stmarm(kflc(1,igc)) ! remove data 00768 n=n+n 00769 IF(nst(1,igc) == 1) THEN 00770 n=n+1 00771 array1(n)=y1 00772 nst(1,igc)=0 ! reset 00773 END IF 00774 IF(n /= 0) THEN 00775 xyplws(7,igc)=xyplws( 9,igc) 00776 xyplws(8,igc)=xyplws(10,igc) 00777 CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc)) 00778 wght=FLOAT(n)/FLOAT(nst(3,igc)*kflc(5,igc)) 00779 xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght 00780 xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght 00781 xyplws(2,igc)=xyplws(5,igc)/(FLOAT(nst(2,igc))+wght) 00782 xyplws(4,igc)=xyplws(6,igc)/(FLOAT(nst(2,igc))+wght) 00783 xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc)) 00784 xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc)) 00785 CALL stmadp(jflc(1,igc),xyplws(1,igc)) 00786 END IF 00787 00788 WRITE(*,*) ' ' 00789 WRITE(*,*) 'Store ',igc,': ',gtext(igc) 00790 IF(jflc(4,igc) == 0) THEN 00791 WRITE(*,*) ' stored n-tuples: ',jflc(3,igc) 00792 ELSE 00793 WRITE(*,*) ' stored n-tuples, not-stored n-tuples: ', & 00794 jflc(3,igc),', ',jflc(4,igc) 00795 END IF 00796 00797 CALL stmacp(jflc(1,igc),array,na) ! get all data 00798 na=na/2 00799 DO n=1,na 00800 WRITE(*,102) n,(array4(j,n),j=1,4) 00801 END DO 00802 END IF 00803 END DO 00804 RETURN 00805 00806 ENTRY gmplun(lunw) ! unit for output 00807 lun=lunw 00808 RETURN 00809 00810 ENTRY gmpwrt(ig) ! write XY text file 00811 IF(lun <= 0) RETURN 00812 IF(ig == 0) THEN 00813 iga=1 00814 igb=numgxy 00815 ELSE 00816 IF(ig <= 0.OR.ig > numgxy) RETURN 00817 iga=ig 00818 igb=ig 00819 END IF 00820 DO igc=iga,igb 00821 IF(igtp(igc) == 5) THEN 00822 00823 CALL stmacp(kflc(1,igc),array,n) ! get data 00824 CALL stmarm(kflc(1,igc)) ! remove data 00825 n=n+n 00826 IF(nst(1,igc) == 1) THEN 00827 n=n+1 00828 array1(n)=y1 00829 nst(1,igc)=0 ! reset 00830 END IF 00831 IF(n /= 0) THEN 00832 xyplws(7,igc)=xyplws( 9,igc) 00833 xyplws(8,igc)=xyplws(10,igc) 00834 CALL rmesig(array1,n,xyplws(2,igc),xyplws(4,igc)) 00835 wght=FLOAT(n)/FLOAT(nst(3,igc)*kflc(5,igc)) 00836 xyplws(5,igc)=xyplws(5,igc)+xyplws(2,igc)*wght 00837 xyplws(6,igc)=xyplws(6,igc)+xyplws(4,igc)*wght 00838 xyplws(2,igc)=xyplws(5,igc)/(FLOAT(nst(2,igc))+wght) 00839 xyplws(4,igc)=xyplws(6,igc)/(FLOAT(nst(2,igc))+wght) 00840 xyplws(1,igc)=0.5*(xyplws(7,igc)+xyplws(8,igc)) 00841 xyplws(3,igc)=0.5*(xyplws(8,igc)-xyplws(7,igc)) 00842 CALL stmadp(jflc(1,igc),xyplws(1,igc)) 00843 END IF 00844 00845 END IF 00846 IF(jflc(3,igc)+jflc(4,igc) /= 0) THEN 00847 WRITE(lun,204) ' ' 00848 WRITE(lun,201) igc,lvers(igc),igtp(igc) 00849 WRITE(lun,204) gtext(igc) 00850 WRITE(lun,203) jflc(3,igc)+jflc(4,igc) 00851 CALL stmacp(jflc(1,igc),array,na) ! get all data 00852 IF(igtp(igc) >= 1.AND.igtp(igc) <= 3) THEN 00853 WRITE(lun,204) 'x-y' 00854 DO n=1,na 00855 WRITE(lun,205) array(1,n),array(2,n) 00856 END DO 00857 ELSE IF(igtp(igc) == 4.OR.igtp(igc) == 5) THEN 00858 WRITE(lun,204) 'x-y-dx-dy' 00859 na=na/2 00860 DO n=1,na 00861 WRITE(lun,205) (array4(j,n),j=1,4) 00862 END DO 00863 END IF 00864 WRITE(lun,204) 'end of xy-data' 00865 END IF 00866 END DO 00867 102 FORMAT(i12,4G15.7) 00868 ! 103 FORMAT(' Index ___X___ ___Y___ '/ 00869 ! + ' ----- -------------- --------------') 00870 ! 104 FORMAT(' Index ___X___ ___Y___ ', 00871 ! + ' ___DX__ ___DY__ '/ 00872 ! + ' ----- -------------- --------------', 00873 ! + ' -------------- --------------') 00874 201 FORMAT('XY-Data ',i4,10X,'version ',i4,10X,'type',i2) 00875 203 FORMAT(10X,'stored not-stored ',2I10) 00876 204 FORMAT(a) 00877 205 FORMAT(3X,4G15.7) 00878 END SUBROUTINE gmpdef 00879 00880 00881 SUBROUTINE stmars ! init/reset storage 00882 IMPLICIT NONE 00883 INTEGER:: i 00884 INTEGER:: ifre 00885 INTEGER:: ifrea 00886 INTEGER:: ifreb 00887 INTEGER:: ind 00888 INTEGER:: j 00889 INTEGER:: n 00890 INTEGER:: ndim 00891 REAL:: x 00892 REAL:: y 00893 PARAMETER (ndim=5000) ! storage dimension, should be NUMGXY*NLIMIT 00894 REAL:: tk(2,ndim) ! pair storage for data pairs 00895 INTEGER:: next(ndim) ! pointer 00896 INTEGER:: iflc1 ! first and last index of free pairs 00897 INTEGER::iflc2 ! first and last index of free pairs 00898 SAVE 00899 00900 REAL:: four(4) ! double_pair, copy array 00901 REAL::array(2,*) ! double_pair, copy array 00902 INTEGER:: jflc(5) ! user array 00903 ! JFLC(1) = first used index 00904 ! JFLC(2) = last used index 00905 ! JFLC(3) = counter of used places 00906 ! JFLC(4) = counter of ignored 00907 ! JFLC(5) = limit for JFLC(3) 00908 ! ... 00909 DO i=1,ndim 00910 next(i)=i+1 ! pointer to next free location 00911 tk(1,i)=0.0 ! reset 00912 tk(2,i)=0.0 00913 END DO 00914 next(ndim)=0 ! ... and end pointer 00915 iflc1=1 ! index first free pair 00916 iflc2=ndim ! index last free pair 00917 RETURN 00918 00919 ENTRY stmapr(jflc,x,y) ! store pair (X,Y) 00920 ifre=iflc1 ! index of free place 00921 IF(ifre == 0.OR.jflc(3) >= jflc(5)) THEN ! overflow 00922 jflc(4)=jflc(4)+1 00923 ELSE 00924 iflc1=next(ifre) ! pointer to new free location 00925 IF(jflc(1) == 0) THEN ! first item 00926 jflc(1)=ifre 00927 ELSE 00928 next(jflc(2))=ifre 00929 END IF 00930 next(ifre)=0 00931 jflc(2)=ifre ! last index 00932 jflc(3)=jflc(3)+1 ! counter 00933 tk(1,ifre)=x 00934 tk(2,ifre)=y 00935 END IF 00936 RETURN 00937 00938 ENTRY stmadp(jflc,four) ! store double pair 00939 ifrea=iflc1 ! index of 1. free place 00940 IF(ifrea == 0) THEN ! overflow 00941 jflc(4)=jflc(4)+1 00942 ELSE 00943 ifreb=next(iflc1) ! index of 2. free place 00944 IF(ifreb == 0.OR.jflc(3) >= 2*jflc(5)) THEN ! overflow 00945 jflc(4)=jflc(4)+1 00946 ELSE 00947 iflc1=next(ifreb) ! pointer to new free location 00948 IF(jflc(1) == 0) THEN ! first item 00949 jflc(1)=ifrea 00950 ELSE 00951 next(jflc(2))=ifrea 00952 END IF 00953 next(ifreb)=0 00954 jflc(2)=ifreb ! last index 00955 jflc(3)=jflc(3)+1 ! counter 00956 tk(1,ifrea)=four(1) 00957 tk(2,ifrea)=four(2) 00958 tk(1,ifreb)=four(3) 00959 tk(2,ifreb)=four(4) 00960 END IF 00961 END IF 00962 RETURN 00963 00964 ENTRY stmacp(jflc,array,n) ! copy (cp) all pairs to array 00965 n=0 00966 ind=jflc(1) 00967 10 IF(ind == 0) RETURN 00968 n=n+1 00969 array(1,n)=tk(1,ind) 00970 array(2,n)=tk(2,ind) 00971 ind=next(ind) 00972 GO TO 10 00973 00974 ENTRY stmarm(jflc) ! remove (rm) stored paiirs 00975 next(iflc2)=jflc(1) ! connect to free space 00976 iflc2=jflc(2) ! new last free index 00977 DO j=1,4 00978 jflc(j)=0 00979 END DO 00980 END SUBROUTINE stmars ! init/ 00981 00982 SUBROUTINE rmesig(x,n,xloc,xsca) ! robust mean and sigma 00983 IMPLICIT NONE 00984 INTEGER:: i 00985 ! robust determination of location and scale parameter, 00986 ! for Gaussian data: location=mean and scale=standard deviation 00987 ! XLOC = median of X_i (N values in array X(N)) 00988 ! XCSA = median of | X_i - XLOC |, times 1.4826 00989 00990 REAL, INTENT(IN OUT) :: x(n) ! input array, modified 00991 INTEGER, INTENT(IN) :: n 00992 REAL, INTENT(OUT) :: xloc 00993 REAL, INTENT(OUT) :: xsca 00994 SAVE 00995 ! ... 00996 xloc=0.0 00997 xsca=0.0 00998 IF(n <= 0) RETURN 00999 CALL heapf(x,n) ! sort 01000 xloc=0.5*(x((n+1)/2)+x((n+2)/2)) ! location 01001 DO i=1,n 01002 x(i)=ABS(x(i)-xloc) 01003 END DO 01004 CALL heapf(x,n) ! sort 01005 xsca=1.4826*0.5*(x((n+1)/2)+x((n+2)/2)) ! dispersion 01006 END SUBROUTINE rmesig 01007 01008 01009
1.8.0
