-
Notifications
You must be signed in to change notification settings - Fork 4
/
mncntr.F
158 lines (158 loc) · 5.07 KB
/
mncntr.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
*
* $Id: mncntr.F,v 1.1.1.1 2000/06/08 11:19:19 andras Exp $
*
* $Log: mncntr.F,v $
* Revision 1.1.1.1 2000/06/08 11:19:19 andras
* import of MINUIT from CERNlib 2000
*
* Revision 1.1.1.1 1996/03/07 14:31:28 mclareni
* Minuit
*
*
#include "minuit/pilot.h"
SUBROUTINE MNCNTR(FCN,KE1,KE2,IERRF,FUTIL)
#include "minuit/d506dp.inc"
CC to print function contours in two variables, on line printer
CC
#include "minuit/d506cm.inc"
EXTERNAL FCN,FUTIL
PARAMETER (NUMBCS=20,NXMAX=115)
DIMENSION CONTUR(NUMBCS), FCNA(NXMAX),FCNB(NXMAX)
CHARACTER CLABEL*(NUMBCS)
CHARACTER CHLN*(NXMAX),CHMID*(NXMAX),CHZERO*(NXMAX)
DATA CLABEL/'0123456789ABCDEFGHIJ'/
C input arguments: parx, pary, devs, ngrid
IF (KE1.LE.0 .OR. KE2.LE.0) GO TO 1350
IF (KE1.GT.NU .OR. KE2.GT.NU) GO TO 1350
KI1 = NIOFEX(KE1)
KI2 = NIOFEX(KE2)
IF (KI1.LE.0 .OR. KI2.LE.0) GO TO 1350
IF (KI1 .EQ. KI2) GO TO 1350
C
IF (ISW(2) .LT. 1) THEN
CALL MNHESS(FCN,FUTIL)
CALL MNWERR
ENDIF
NPARX = NPAR
XSAV = U(KE1)
YSAV = U(KE2)
DEVS = WORD7(3)
IF (DEVS .LE. ZERO) DEVS=2.
XLO = U(KE1) - DEVS*WERR(KI1)
XUP = U(KE1) + DEVS*WERR(KI1)
YLO = U(KE2) - DEVS*WERR(KI2)
YUP = U(KE2) + DEVS*WERR(KI2)
NGRID = WORD7(4)
IF (NGRID .LE. 0) THEN
NGRID=25
NX = MIN(NPAGWD-15,NGRID)
NY = MIN(NPAGLN-7, NGRID)
ELSE
NX = NGRID
NY = NGRID
ENDIF
IF (NX .LT. 11) NX=11
IF (NY .LT. 11) NY=11
IF (NX .GE. NXMAX) NX=NXMAX-1
C ask if parameter outside limits
IF (NVARL(KE1) .GT. 1) THEN
IF (XLO .LT. ALIM(KE1)) XLO = ALIM(KE1)
IF (XUP .GT. BLIM(KE1)) XUP = BLIM(KE1)
ENDIF
IF (NVARL(KE2) .GT. 1) THEN
IF (YLO .LT. ALIM(KE2)) YLO = ALIM(KE2)
IF (YUP .GT. BLIM(KE2)) YUP = BLIM(KE2)
ENDIF
BWIDX = (XUP-XLO)/REAL(NX)
BWIDY = (YUP-YLO)/REAL(NY)
IXMID = INT((XSAV-XLO)*REAL(NX)/(XUP-XLO)) + 1
IF (AMIN .EQ. UNDEFI) CALL MNAMIN(FCN,FUTIL)
DO 185 I= 1, NUMBCS
CONTUR(I) = AMIN + UP*FLOAT(I-1)**2
185 CONTINUE
CONTUR(1) = CONTUR(1) + 0.01*UP
C fill FCNB to prepare first row, and find column zero
U(KE2) = YUP
IXZERO = 0
XB4 = ONE
DO 200 IX= 1, NX+1
U(KE1) = XLO + REAL(IX-1)*BWIDX
CALL FCN(NPARX,GIN,FF,U,4,FUTIL)
FCNB(IX) = FF
IF (XB4.LT.ZERO .AND. U(KE1).GT.ZERO) IXZERO = IX-1
XB4 = U(KE1)
CHMID(IX:IX) = '*'
CHZERO(IX:IX)= '-'
200 CONTINUE
WRITE (ISYSWR,'(A,I3,A,A)') ' Y-AXIS: PARAMETER ',
+ KE2,': ',CPNAM(KE2)
IF (IXZERO .GT. 0) THEN
CHZERO(IXZERO:IXZERO) = '+'
CHLN = ' '
WRITE (ISYSWR,'(12X,A,A)') CHLN(1:IXZERO),'X=0'
ENDIF
C loop over rows
DO 280 IY= 1, NY
UNEXT = U(KE2) - BWIDY
C prepare this line's background pattern for contour
CHLN = ' '
CHLN(IXMID:IXMID) = '*'
IF (IXZERO .NE. 0) CHLN(IXZERO:IXZERO) = ':'
IF (U(KE2).GT.YSAV .AND. UNEXT.LT.YSAV) CHLN=CHMID
IF (U(KE2).GT.ZERO .AND. UNEXT.LT.ZERO) CHLN=CHZERO
U(KE2) = UNEXT
YLABEL = U(KE2) + 0.5*BWIDY
C move FCNB to FCNA and fill FCNB with next row
DO 220 IX= 1, NX+1
FCNA(IX) = FCNB(IX)
U(KE1) = XLO + REAL(IX-1)*BWIDX
CALL FCN(NPARX,GIN,FF,U,4,FUTIL)
FCNB(IX) = FF
220 CONTINUE
C look for contours crossing the FCNxy squares
DO 250 IX= 1, NX
FMX = MAX(FCNA(IX),FCNB(IX),FCNA(IX+1),FCNB(IX+1))
FMN = MIN(FCNA(IX),FCNB(IX),FCNA(IX+1),FCNB(IX+1))
DO 230 ICS= 1, NUMBCS
IF (CONTUR(ICS) .GT. FMN) GO TO 240
230 CONTINUE
GO TO 250
240 IF (CONTUR(ICS) .LT. FMX) CHLN(IX:IX)=CLABEL(ICS:ICS)
250 CONTINUE
C print a row of the contour plot
WRITE (ISYSWR,'(1X,G12.4,1X,A)') YLABEL,CHLN(1:NX)
280 CONTINUE
C contours printed, label x-axis
CHLN = ' '
CHLN( 1: 1) = 'I'
CHLN(IXMID:IXMID) = 'I'
CHLN(NX:NX) = 'I'
WRITE (ISYSWR,'(14X,A)') CHLN(1:NX)
C the hardest of all: print x-axis scale!
CHLN = ' '
IF (NX .LE. 26) THEN
NL = MAX(NX-12,2)
NL2 = NL/2
WRITE (ISYSWR,'(8X,G12.4,A,G12.4)') XLO,CHLN(1:NL),XUP
WRITE (ISYSWR,'(14X,A,G12.4)') CHLN(1:NL2),XSAV
ELSE
NL = MAX(NX-24,2)/2
NL2 = NL
IF (NL .GT. 10) NL2=NL-6
WRITE (ISYSWR,'(8X,G12.4,A,G12.4,A,G12.4)') XLO,
+ CHLN(1:NL),XSAV,CHLN(1:NL2),XUP
ENDIF
WRITE (ISYSWR,'(6X,A,I3,A,A,A,G12.4)') ' X-AXIS: PARAMETER',
+ KE1,': ',CPNAM(KE1),' ONE COLUMN=',BWIDX
WRITE (ISYSWR,'(A,G12.4,A,G12.4,A)') ' FUNCTION VALUES: F(I)=',
+ AMIN,' +',UP,' *I**2'
C finished. reset input values
U(KE1) = XSAV
U(KE2) = YSAV
IERRF = 0
RETURN
1350 WRITE (ISYSWR,1351)
1351 FORMAT (' INVALID PARAMETER NUMBER(S) REQUESTED. IGNORED.' /)
IERRF = 1
RETURN
END