source: trunk/MagicSoft/Simulation/Corsika/Mmcs/amoeba.f

Last change on this file was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 3.7 KB
Line 
1 SUBROUTINE AMOEBA(P,Y,MP,NP,NDIM,FTOL,FUNK,ITER,IFLAG)
2
3C-----------------------------------------------------------------------
4C
5C FITTING ROUTINE
6C REFERENCE : NUMERICAL RECIPES, W.H. PRESS ET AL.,
7C CAMBRIDGE UNIVERSITY PRESS, 1992 ISBN 0 521 43064 X
8C ADAPTED FOR DOUBLE PRECISION
9C THIS SUBROUTINE IS CALLED FROM LONGFT
10C-----------------------------------------------------------------------
11
12 IMPLICIT NONE
13*KEEP,RUNPAR.
14 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
15 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
16 * MONIOU,MDEBUG,NUCNUC,
17 * CETAPE,
18 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
19 * N1STTR,MDBASE,
20 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
21 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
22 * ,GHEISH,GHESIG
23 COMMON /RUNPAC/ DSN,HOST,USER
24 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
25 REAL STEPFC
26 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
27 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
28 * N1STTR,MDBASE
29 INTEGER CETAPE
30 CHARACTER*79 DSN
31 CHARACTER*20 HOST,USER
32
33 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
34 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
35 * ,GHEISH,GHESIG
36*KEND.
37
38 INTEGER ITMAX,MP,NMAX,NP
39C MAXIMUM NUMBER OF TRIAL PER CALL
40 PARAMETER (ITMAX=5000)
41 PARAMETER (NMAX=20)
42 DOUBLE PRECISION AMOTRY,FTOL,FUNK,P(MP,NP),PSUM(NMAX),
43 * RTOL,SUM,SWAP,Y(MP),YSAVE,YTRY
44 INTEGER I,IFLAG,IHI,ILO,INHI,ITER,J,M,N,NDIM
45 EXTERNAL FUNK
46
47CU USES AMOTRY,FUNK
48C-----------------------------------------------------------------------
49
50 IF ( DEBUG ) WRITE(MDEBUG,*) 'AMOEBA:'
51
52 IFLAG = 0
53 ITER = 0
54 1 DO 12 N=1,NDIM
55 SUM = 0.D0
56 DO 11 M=1,NDIM+1
57 SUM = SUM + P(M,N)
58 11 CONTINUE
59 PSUM(N) = SUM
60 12 CONTINUE
61 2 ILO=1
62 IF ( Y(1) .GT. Y(2) ) THEN
63 IHI = 1
64 INHI = 2
65 ELSE
66 IHI = 2
67 INHI = 1
68 ENDIF
69 DO 13 I=1,NDIM+1
70 IF ( Y(I) .LE. Y(ILO) ) ILO = I
71 IF ( Y(I) .GT. Y(IHI) ) THEN
72 INHI = IHI
73 IHI = I
74 ELSEIF ( Y(I) .GT. Y(INHI) ) THEN
75 IF ( I .NE. IHI ) INHI = I
76 ENDIF
77 13 CONTINUE
78 RTOL = 2.D0*ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))
79 IF ( RTOL .LT. FTOL ) THEN
80 SWAP = Y(1)
81 Y(1) = Y(ILO)
82 Y(ILO) = SWAP
83 DO 14 N=1,NDIM
84 SWAP = P(1,N)
85 P(1,N) = P(ILO,N)
86 P(ILO,N) = SWAP
87 14 CONTINUE
88 RETURN
89 ENDIF
90 IF ( ITER .GE.ITMAX ) THEN
91 IF(DEBUG) WRITE(MDEBUG,*) 'AMOEBA: ITMAX EXCEEDED IN AMOEBA'
92 IFLAG = 1
93 RETURN
94 ENDIF
95 ITER = ITER + 2
96 YTRY = AMOTRY(P,Y,PSUM,MP,NP,NDIM,FUNK,IHI,-1.0D0)
97 IF ( YTRY .LE. Y(ILO) ) THEN
98 YTRY = AMOTRY(P,Y,PSUM,MP,NP,NDIM,FUNK,IHI,2.0D0)
99 ELSEIF ( YTRY .GE. Y(INHI) ) THEN
100 YSAVE = Y(IHI)
101 YTRY = AMOTRY(P,Y,PSUM,MP,NP,NDIM,FUNK,IHI,0.5D0)
102 IF ( YTRY .GE. YSAVE ) THEN
103 DO 16 I=1,NDIM+1
104 IF ( I .NE. ILO ) THEN
105 DO 15 J=1,NDIM
106 PSUM(J) = 0.5D0 * (P(I,J) + P(ILO,J))
107 P(I,J) = PSUM(J)
108 15 CONTINUE
109 Y(I) = FUNK(PSUM)
110 ENDIF
111 16 CONTINUE
112 ITER = ITER + NDIM
113 GOTO 1
114 ENDIF
115 ELSE
116 ITER = ITER - 1
117 ENDIF
118 GOTO 2
119 END
120C=======================================================================
121
Note: See TracBrowser for help on using the repository browser.