COSMOS v7.655  COSMOSv7655
(AirShowerMC)
Functions/Subroutines
mncros.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine mncros (FCN, AOPT, IERCR, FUTIL)
 

Function/Subroutine Documentation

◆ mncros()

subroutine mncros ( external  FCN,
  AOPT,
  IERCR,
external  FUTIL 
)

Definition at line 10 of file mncros.f.

References a, b, c, d, e, f, false, g, h, softenpik::half, i, j, m, maxp, mneval(), mnpfit(), mnplot(), mnwarn(), n, o, p, parameter(), true, up(), x, x1(), x2(), xs, and z.

Referenced by mncont(), and mnmnot().

10 *
11 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
12 *
13 * $Log: d506dp.inc,v $
14 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
15 * Minuit
16 *
17 *
18 *
19 *
20 * d506dp.inc
21 *
22 C ************ DOUBLE PRECISION VERSION *************
23  IMPLICIT DOUBLE PRECISION (a-h,o-z)
24 CC Find point where MNEVAL=AMIN+UP, along the line through
25 CC XMIDCR,YMIDCR with direction XDIRCR,YDIRCR, where X and Y
26 CC are parameters KE1CR and KE2CR. If KE2CR=0 (from MINOS),
27 CC only KE1CR is varied. From MNCONT, both are varied.
28 CC Crossing point is at
29 CC (U(KE1),U(KE2)) = (XMID,YMID) + AOPT*(XDIR,YDIR)
30 CC
31 *
32 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
33 *
34 * $Log: d506cm.inc,v $
35 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
36 * Minuit
37 *
38 *
39 *
40 *
41 * d506cm.inc
42 *
43  parameter(mne=100 , mni=50)
44  parameter(mnihl=mni*(mni+1)/2)
45  CHARACTER*10 cpnam
46  COMMON
47  1/mn7nam/ cpnam(mne)
48  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
49  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
50  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
51  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
52  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
53  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
54  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
55  9/mn7fx1/ ipfix(mni) ,npfix
56  a/mn7var/ vhmat(mnihl)
57  b/mn7vat/ vthmat(mnihl)
58  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
59 C
60  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
61  parameter(zero=0.0, one=1.0, half=0.5)
62  COMMON
63  d/mn7npr/ maxint ,npar ,maxext ,nu
64  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
65  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
66  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
67  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
68  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
69  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
70  j/mn7arg/ word7(maxp)
71  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
72  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
73  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
74  n/mn7cpt/ chpt(maxcpt)
75  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
76  CHARACTER ctitl*50, cword*(maxcwd), cundef*10, cfrom*8,
77  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
78  LOGICAL lwarn, lrepor, limset, lnolim, lnewmn, lphead
79  CHARACTER chere*10, charal*28, chsign*4
80  parameter(chere='MNCROS ', mlsb=3, maxitr=15, tlr=0.01)
81  dimension flsb(mlsb),alsb(mlsb), coeff(3)
82  LOGICAL ldebug
83  EXTERNAL fcn,futil
84  DATA charal/' .ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
85  ldebug = (idbg(6) .GE. 1)
86  aminsv = amin
87 C convergence when F is within TLF of AIM and next prediction
88 C of AOPT is within TLA of previous value of AOPT
89  aim = amin + up
90  tlf = tlr*up
91  tla = tlr
92  xpt(1) = 0.0
93  ypt(1) = aim
94  chpt(1) = ' '
95  ipt = 1
96  IF (ke2cr .EQ. 0) THEN
97  xpt(2) = -1.0
98  ypt(2) = amin
99  chpt(2) = '.'
100  ipt = 2
101  ENDIF
102 C find the largest allowed A
103  aulim = 100.
104  DO 100 ik= 1, 2
105  IF (ik .EQ. 1) THEN
106  kex = ke1cr
107  zmid = xmidcr
108  zdir = xdircr
109  ELSE
110  IF (ke2cr .EQ. 0) GO TO 100
111  kex = ke2cr
112  zmid = ymidcr
113  zdir = ydircr
114  ENDIF
115  IF (nvarl(kex) .LE. 1) GO TO 100
116  IF (zdir .EQ. zero) GO TO 100
117  zlim = alim(kex)
118  IF (zdir .GT. zero) zlim = blim(kex)
119  aulim = min(aulim,(zlim-zmid)/zdir)
120  100 CONTINUE
121 C LSB = Line Search Buffer
122 C first point
123  anext = 0.
124  aopt = anext
125  limset = .false.
126  IF (aulim .LT. aopt+tla) limset = .true.
127  CALL mneval(fcn,anext,fnext,ierev,futil)
128 C debug printout:
129  IF (ldebug) WRITE (isyswr,'(A,I8,A,F10.5,A,2F10.5)')
130  + ' MNCROS: calls=',nfcn,' AIM=',aim,' F,A=',fnext,aopt
131  IF (ierev .GT. 0) GO TO 900
132  IF (limset .AND. fnext .LE. aim) GO TO 930
133  ipt = ipt + 1
134  xpt(ipt) = anext
135  ypt(ipt) = fnext
136  chpt(ipt)= charal(ipt:ipt)
137  alsb(1) = anext
138  flsb(1) = fnext
139  fnext = max(fnext,aminsv+0.1*up)
140  aopt = sqrt((up)/(fnext-aminsv)) - 1.0
141  IF (abs(fnext-aim) .LT. tlf) GO TO 800
142 C
143  IF (aopt .LT. -half) aopt = -half
144  IF (aopt .GT. one) aopt = one
145  limset = .false.
146  IF (aopt .GT. aulim) THEN
147  aopt = aulim
148  limset = .true.
149  ENDIF
150  CALL mneval(fcn,aopt,fnext,ierev,futil)
151 C debug printout:
152  IF (ldebug) WRITE (isyswr,'(A,I8,A,F10.5,A,2F10.5)')
153  + ' MNCROS: calls=',nfcn,' AIM=',aim,' F,A=',fnext,aopt
154  IF (ierev .GT. 0) GO TO 900
155  IF (limset .AND. fnext .LE. aim) GO TO 930
156  alsb(2) = aopt
157  ipt = ipt + 1
158  xpt(ipt) = alsb(2)
159  ypt(ipt) = fnext
160  chpt(ipt)= charal(ipt:ipt)
161  flsb(2) = fnext
162  dfda = (flsb(2)-flsb(1))/ (alsb(2)-alsb(1))
163 C DFDA must be positive on the contour
164  IF (dfda .GT. zero) GO TO 460
165  300 CALL mnwarn('D',chere,'Looking for slope of the right sign')
166  maxlk = maxitr - ipt
167  DO 400 it= 1, maxlk
168  alsb(1) = alsb(2)
169  flsb(1) = flsb(2)
170  aopt = alsb(1) + 0.2*REAL(it)
171  limset = .false.
172  IF (aopt .GT. aulim) THEN
173  aopt = aulim
174  limset = .true.
175  ENDIF
176  CALL mneval(fcn,aopt,fnext,ierev,futil)
177 C debug printout:
178  IF (ldebug) WRITE (isyswr,'(A,I8,A,F10.5,A,2F10.5)')
179  + ' MNCROS: calls=',nfcn,' AIM=',aim,' F,A=',fnext,aopt
180  IF (ierev .GT. 0) GO TO 900
181  IF (limset .AND. fnext .LE. aim) GO TO 930
182  alsb(2) = aopt
183  ipt = ipt + 1
184  xpt(ipt) = alsb(2)
185  ypt(ipt) = fnext
186  chpt(ipt)= charal(ipt:ipt)
187  flsb(2) = fnext
188  dfda = (flsb(2)-flsb(1))/ (alsb(2)-alsb(1))
189  IF (dfda .GT. zero) GO TO 450
190  400 CONTINUE
191  CALL mnwarn('W',chere,'Cannot find slope of the right sign')
192  GO TO 950
193  450 CONTINUE
194 C we have two points with the right slope
195  460 aopt = alsb(2) + (aim-flsb(2))/dfda
196  fdist = min(abs(aim -flsb(1)),abs(aim -flsb(2)))
197  adist = min(abs(aopt-alsb(1)),abs(aopt-alsb(2)))
198  tla = tlr
199  IF (abs(aopt) .GT. one) tla = tlr*abs(aopt)
200  IF (adist .LT. tla .AND. fdist .LT. tlf) GO TO 800
201  IF (ipt .GE. maxitr) GO TO 950
202  bmin = min(alsb(1),alsb(2)) - 1.0
203  IF (aopt .LT. bmin) aopt = bmin
204  bmax = max(alsb(1),alsb(2)) + 1.0
205  IF (aopt .GT. bmax) aopt = bmax
206 C Try a third point
207  limset = .false.
208  IF (aopt .GT. aulim) THEN
209  aopt = aulim
210  limset = .true.
211  ENDIF
212  CALL mneval(fcn,aopt,fnext,ierev,futil)
213 C debug printout:
214  IF (ldebug) WRITE (isyswr,'(A,I8,A,F10.5,A,2F10.5)')
215  + ' MNCROS: calls=',nfcn,' AIM=',aim,' F,A=',fnext,aopt
216  IF (ierev .GT. 0) GO TO 900
217  IF (limset .AND. fnext .LE. aim) GO TO 930
218  alsb(3) = aopt
219  ipt = ipt + 1
220  xpt(ipt) = alsb(3)
221  ypt(ipt) = fnext
222  chpt(ipt)= charal(ipt:ipt)
223  flsb(3) = fnext
224  inew = 3
225 C now we have three points, ask how many <AIM
226  ecarmn = abs(fnext-aim)
227  ibest = 3
228  ecarmx = 0.
229  noless = 0
230  DO 480 i= 1, 3
231  ecart = abs(flsb(i) - aim)
232  IF (ecart .GT. ecarmx) THEN
233  ecarmx = ecart
234  iworst = i
235  ENDIF
236  IF (ecart .LT. ecarmn) THEN
237  ecarmn = ecart
238  ibest = i
239  ENDIF
240  IF (flsb(i) .LT. aim) noless = noless + 1
241  480 CONTINUE
242  inew = ibest
243 C if at least one on each side of AIM, fit a parabola
244  IF (noless.EQ.1 .OR. noless.EQ.2) GO TO 500
245 C if all three are above AIM, third must be closest to AIM
246  IF (noless .EQ. 0 .AND. ibest .NE. 3) GO TO 950
247 C if all three below, and third is not best, then slope
248 C has again gone negative, look for positive slope.
249  IF (noless .EQ. 3 .AND. ibest .NE. 3) THEN
250  alsb(2) = alsb(3)
251  flsb(2) = flsb(3)
252  GO TO 300
253  ENDIF
254 C in other cases, new straight line thru last two points
255  alsb(iworst) = alsb(3)
256  flsb(iworst) = flsb(3)
257  dfda = (flsb(2)-flsb(1))/ (alsb(2)-alsb(1))
258  GO TO 460
259 C parabola fit
260  500 CALL mnpfit(alsb,flsb,3,coeff,sdev)
261  IF (coeff(3) .LE. zero) CALL mnwarn ('D',chere,
262  + 'Curvature is negative near contour line.')
263  determ = coeff(2)**2 - 4.*coeff(3)*(coeff(1)-aim)
264  IF (determ .LE. zero) THEN
265  CALL mnwarn('D',chere,'Problem 2, impossible determinant')
266  GO TO 950
267  ENDIF
268 C Find which root is the right one
269  rt = sqrt(determ)
270  x1 = (-coeff(2) + rt)/(2.*coeff(3))
271  x2 = (-coeff(2) - rt)/(2.*coeff(3))
272  s1 = coeff(2) + 2.*x1*coeff(3)
273  s2 = coeff(2) + 2.*x2*coeff(3)
274  IF (s1*s2 .GT. zero) WRITE (isyswr,'(A)') ' MNCONTour problem 1'
275  aopt = x1
276  slope = s1
277  IF (s2 .GT. zero) THEN
278  aopt = x2
279  slope = s2
280  ENDIF
281 C ask if converged
282  tla = tlr
283  IF (abs(aopt) .GT. one) tla = tlr*abs(aopt)
284  IF (abs(aopt-alsb(ibest)) .LT. tla .AND.
285  & abs(flsb(ibest)-aim) .LT. tlf) GO TO 800
286  IF (ipt .GE. maxitr) GO TO 950
287 C see if proposed point is in acceptable zone between L and R
288 C first find ILEFT, IRIGHT, IOUT and IBEST
289  ileft = 0
290  iright = 0
291  ibest = 1
292  ecarmx = 0.
293  ecarmn = abs(aim-flsb(1))
294  DO 550 i= 1, 3
295  ecart = abs(flsb(i) - aim)
296  IF (ecart .LT. ecarmn) THEN
297  ecarmn = ecart
298  ibest = i
299  ENDIF
300  IF (ecart .GT. ecarmx) ecarmx = ecart
301  IF (flsb(i) .GT. aim) THEN
302  IF (iright .EQ. 0) THEN
303  iright = i
304  ELSE IF (flsb(i) .GT. flsb(iright)) THEN
305  iout = i
306  ELSE
307  iout = iright
308  iright = i
309  ENDIF
310  ELSE IF (ileft .EQ. 0) THEN
311  ileft = i
312  ELSE IF (flsb(i) .LT. flsb(ileft)) THEN
313  iout = i
314  ELSE
315  iout = ileft
316  ileft = i
317  ENDIF
318  550 CONTINUE
319 C avoid keeping a very bad point next time around
320  IF (ecarmx .GT. 10.*abs(flsb(iout)-aim))
321  & aopt = half*aopt + half*half*(alsb(iright)+alsb(ileft))
322 C knowing ILEFT and IRIGHT, get acceptable window
323  smalla = 0.1*tla
324  IF (slope*smalla .GT. tlf) smalla = tlf/slope
325  aleft = alsb(ileft) + smalla
326  aright = alsb(iright) - smalla
327 C move proposed point AOPT into window if necessary
328  IF (aopt .LT. aleft) aopt = aleft
329  IF (aopt .GT. aright) aopt = aright
330  IF (aleft .GT. aright)aopt = half*(aleft + aright)
331 C see if proposed point outside limits (should be impossible!)
332  limset = .false.
333  IF (aopt .GT. aulim) THEN
334  aopt = aulim
335  limset = .true.
336  ENDIF
337 C Evaluate function at new point AOPT
338  CALL mneval(fcn,aopt,fnext,ierev,futil)
339 C debug printout:
340  IF (ldebug) WRITE (isyswr,'(A,I8,A,F10.5,A,2F10.5)')
341  + ' MNCROS: calls=',nfcn,' AIM=',aim,' F,A=',fnext,aopt
342  IF (ierev .GT. 0) GO TO 900
343  IF (limset .AND. fnext .LE. aim) GO TO 930
344  ipt = ipt + 1
345  xpt(ipt) = aopt
346  ypt(ipt) = fnext
347  chpt(ipt)= charal(ipt:ipt)
348 C Replace odd point by new one
349  alsb(iout) = aopt
350  flsb(iout) = fnext
351 C the new point may not be the best, but it is the only one
352 C which could be good enough to pass convergence criteria
353  ibest = iout
354  GO TO 500
355 C
356 C Contour has been located, return point to MNCONT OR MINOS
357  800 CONTINUE
358  iercr = 0
359  GO TO 1000
360 C error in the minimization
361  900 IF (ierev .EQ. 1) GO TO 940
362  GO TO 950
363 C parameter up against limit
364  930 iercr = 1
365  GO TO 1000
366 C too many calls to FCN
367  940 iercr = 2
368  GO TO 1000
369 C cannot find next point
370  950 iercr = 3
371 C in any case
372  1000 CONTINUE
373  IF (ldebug) THEN
374  itoohi = 0
375  DO 1100 i= 1, ipt
376  IF (ypt(i) .GT. aim+up) THEN
377  ypt(i) = aim+up
378  chpt(i) = '+'
379  itoohi = 1
380  ENDIF
381  1100 CONTINUE
382  chsign = 'POSI'
383  IF (xdircr .LT. zero) chsign = 'NEGA'
384  IF (ke2cr .EQ. 0) WRITE (isyswr, '(2X,A,A,I3)')
385  + chsign,'TIVE MINOS ERROR, PARAMETER ',ke1cr
386  IF (itoohi .EQ. 1) WRITE (isyswr, '(10X,A)')
387  + 'POINTS LABELLED "+" WERE TOO HIGH TO PLOT.'
388  IF (iercr .EQ. 1) WRITE (isyswr,'(10X,A)')
389  + 'RIGHTMOST POINT IS UP AGAINST LIMIT.'
390  CALL mnplot(xpt,ypt,chpt,ipt,isyswr,npagwd,npagln)
391  ENDIF
392  RETURN
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
z
nodes z
Definition: cstdatmos0.fNew2.h:93
x1
block data include Zlatfit h c fitting region data x1(1)/0.03/
up
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data up(2, 1)/7.0d0/
mneval
subroutine mneval(FCN, ANEXT, FNEXT, IEREV, FUTIL)
Definition: mneval.f:10
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
i
nodes i
Definition: cstdatmos0.fNew2.h:93
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
xs
! constants thru Cosmos real ! if multiplied to deg radian Torad ! light velocity m sec ! infinty ! kg m2 *Togpcm2 g cm2 ! g cm2 *Tokgpm2 kg m2 ! cm *Tom m ! m *Tocm cm ! g cm3 *Tokgpm3 kg m3 ! kg m3 *Togpcm3 g cm3 ! sec *Tonsec nsec ! Tesla m ! Avogadro *A2deninv ! mfp *n * xs
Definition: Zglobalc.h:18
p
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
o
struct ob o[NpMax]
Definition: Zprivate.h:34
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
mnpfit
subroutine mnpfit(PARX2P, PARY2P, NPAR2P, COEF2P, SDEV2P)
Definition: mnpfit.f:10
m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
false
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
maxp
!onst int maxp
Definition: Zprivate.h:3
x2
block data include Zlatfit h c fitting region data x2(1)/0.5/data x1(2)/0.3/
mnwarn
subroutine mnwarn(COPT, CORG, CMES)
Definition: mnwarn.f:10
n
integer n
Definition: Zcinippxc.h:1
softenpik::half
integer, parameter half
Definition: csoftenPiK.f:108
mnplot
subroutine mnplot(XPT, YPT, CHPT, NXYPT, NUNIT, NPAGWD, NPAGLN)
Definition: mnplot.f:10
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
f
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function: