c Program bestfit060318 (KSKLum, 3/18/06) c This program finds the best-fit three-parameter function to a set of c data points. implicit none real AAA, BBB integer NPOINTS, NGRIDPOINTS parameter (AAA = 4.375) parameter (BBB = 2.25) parameter (NPOINTS = 17) c 15 x 15 x 15 grid parameter (NGRIDPOINTS = 15) real x(NPOINTS), y(NPOINTS), rmsdevmin, yyybestfit, 1 rmsdev(NGRIDPOINTS, NGRIDPOINTS, NGRIDPOINTS), 1 avalue(NGRIDPOINTS, NGRIDPOINTS, NGRIDPOINTS), 1 bvalue(NGRIDPOINTS, NGRIDPOINTS, NGRIDPOINTS), 1 cvalue(NGRIDPOINTS, NGRIDPOINTS, NGRIDPOINTS), 1 cellsizea, cellsizeb, cellsizec, acenter, bcenter, ccenter, 1 amin, amax, bmin, bmax, cmin, cmax, aspan, bspan, cspan, 1 rmsdevminlast, avaluelast, bvaluelast, cvaluelast integer imin, jmin, kmin, i, j, k, l, iargc, nargs, itrack, 1 ifound, iterations, maxiterations, nsame character argv*100 data x / - 1., - 0.875, - 0.75, - 0.625, - 0.5, - 0.375, - 0.25, 1 - 0.125, 0., 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1. / data y / 2.25, 2.5, 3.5, 3.75, 4., 4., 4.25, 4.375, 4.5, 4.375, 1 4.25, 4., 4., 3.75, 3.5, 2.5, 2.25 / nargs = 7 if (iargc() .eq. nargs) then open (1, file = 'bestfit060318.tmp') do i = 1, nargs call getarg(i, argv) write (1, *) argv enddo close (1) open (1, file = 'bestfit060318.tmp', status = 'old') read (1, *) amin read (1, *) amax read (1, *) bmin read (1, *) bmax read (1, *) cmin read (1, *) cmax read (1, *) maxiterations close (1) else print *, '*** YOU HAVE ENTERED AN INCORRECT COMMAND LINE ***' print *, 'Correct Command Line:' print *, 'bestfit060318.e amin amax bmin bmax cmin cmax' print *, 'maxiterations' print *, 'amin = minimum value for parameter A' print *, 'amax = maximum value for parameter A' print *, 'bmin = minimum value for parameter B' print *, 'bmax = maximum value for parameter B' print *, 'cmin = minimum value for parameter C' print *, 'cmax = maximum value for parameter C' print *, 'maxiterations = maximum number of iterations' call flush(6) call exit(-1) endif c Modify x and y values. do i = 1, NPOINTS x(i) = x(i) * AAA y(i) = y(i) - BBB enddo c Build grid of A, B, and C values. acenter = 0.5 * (amin + amax) bcenter = 0.5 * (bmin + bmax) ccenter = 0.5 * (cmin + cmax) aspan = amax - amin bspan = bmax - bmin cspan = cmax - cmin cellsizea = (amax - amin) / float(NGRIDPOINTS - 1) cellsizeb = (bmax - bmin) / float(NGRIDPOINTS - 1) cellsizec = (cmax - cmin) / float(NGRIDPOINTS - 1) open (1, file = 'bestfit060318.txt') ifound = 0 iterations = 0 rmsdevminlast = 1.e+6 avaluelast = 1.e+6 bvaluelast = 1.e+6 cvaluelast = 1.e+6 nsame = 0 do while ((ifound .eq. 0) .and. (iterations .lt. maxiterations) 1 .and. (nsame .lt. 20)) iterations = iterations + 1 amin = acenter - 0.5 * aspan amax = acenter + 0.5 * aspan bmin = bcenter - 0.5 * bspan bmax = bcenter + 0.5 * bspan cmin = ccenter - 0.5 * cspan cmax = ccenter + 0.5 * cspan do i = 1, NGRIDPOINTS do j = 1, NGRIDPOINTS do k = 1, NGRIDPOINTS if (NGRIDPOINTS .gt. 1) then avalue(i, j, k) = amin + float(i - 1) * (amax - amin) / 1 float(NGRIDPOINTS - 1) bvalue(i, j, k) = bmin + float(j - 1) * (bmax - bmin) / 1 float(NGRIDPOINTS - 1) cvalue(i, j, k) = cmin + float(k - 1) * (cmax - cmin) / 1 float(NGRIDPOINTS - 1) else if (NGRIDPOINTS .eq. 1) then avalue(i, j, k) = 0.5 * (amin + amax) bvalue(i, j, k) = 0.5 * (bmin + bmax) cvalue(i, j, k) = 0.5 * (cmin + cmax) endif enddo enddo enddo c Calculate RMS deviations. do i = 1, NGRIDPOINTS do j = 1, NGRIDPOINTS do k = 1, NGRIDPOINTS rmsdev(i, j, k) = 0. itrack = 0 do l = 1, NPOINTS if (abs(x(l) / avalue(i, j, k)) .lt. 1.) then yyybestfit = bvalue(i, j, k) * 1 sqrt(1. - (x(l) / avalue(i, j, k))**2.) 1 - cvalue(i, j, k) rmsdev(i, j, k) = rmsdev(i, j, k) + 1 (yyybestfit - y(l))**2. itrack = itrack + 1 endif enddo if (itrack .eq. NPOINTS) then rmsdev(i, j, k) = sqrt(rmsdev(i, j, k) / float(itrack)) else rmsdev(i, j, k) = - 1. endif enddo enddo enddo c Find lowest RMS deviation. rmsdevmin = 1.e+6 do i = 1, NGRIDPOINTS write (1, 100) iterations, avalue(i, 1, 1) 100 format ('iterations, avalue = ', i10, f10.4) do j = 1, NGRIDPOINTS do k = 1, NGRIDPOINTS if ((rmsdev(i, j, k) .lt. rmsdevmin) .and. 1 (rmsdev(i, j, k) .ge. 0.)) then rmsdevmin = rmsdev(i, j, k) imin = i jmin = j kmin = k endif enddo write (1, 110) (rmsdev(i, j, k), k = 1, NGRIDPOINTS) 110 format (f10.4) enddo enddo print *, 'iterations, rmsdevmin, a, b, c = ', 1 iterations, rmsdevmin, avalue(imin, jmin, kmin), 1 bvalue(imin, jmin, kmin), cvalue(imin, jmin, kmin) call flush(6) c Stop program if solution has converged. if ((abs(rmsdevmin - rmsdevminlast) .lt. 1.e-12) .and. 1 (abs(avalue(imin, jmin, kmin) - avaluelast) .lt. 1.e-12) .and. 1 (abs(bvalue(imin, jmin, kmin) - bvaluelast) .lt. 1.e-12) .and. 1 (abs(cvalue(imin, jmin, kmin) - cvaluelast) .lt. 1.e-12)) then nsame = nsame + 1 else nsame = 0 endif rmsdevminlast = rmsdevmin avaluelast = avalue(imin, jmin, kmin) bvaluelast = bvalue(imin, jmin, kmin) cvaluelast = cvalue(imin, jmin, kmin) c Modify acenter, bcenter, and ccenter, and reduce aspan, bspan, and c cspan. Add 0.5 * aspan / float(NGRIDPOINTS - 1) etc. to prevent c program from getting stuck at one grid point. acenter = avalue(imin, jmin, kmin) + 1 0.5 * aspan / float(NGRIDPOINTS - 1) bcenter = bvalue(imin, jmin, kmin) + 1 0.5 * bspan / float(NGRIDPOINTS - 1) ccenter = cvalue(imin, jmin, kmin) + 1 0.5 * cspan / float(NGRIDPOINTS - 1) aspan = 0.95 * aspan bspan = 0.95 * bspan cspan = 0.95 * cspan enddo ! end of "do while (ifound .eq. 0)" close (1) stop end