subroutine FCN (npar, grad, chi2, par, iflag, futil) c Author: Glen Cowan c Date: 8-Dec-1998 c Generic FCN routine for use with MINUIT c Input: integer npar number of parameters to fit c double precision par(npar) parameter vector c integer iflag select what to do c double precision futil optional external function c c Output: double precision grad(npar) gradient vector (not filled) c double precision chi2 chi^2 to be minimized c Use KUIP routines to read in values: c KUPROR for real, KUPROI for integer, KUPROS for character string implicit NONE c constants integer max_points parameter (max_points = 1000) integer max_sets parameter (max_sets = 100) integer nf parameter (nf = 180) ! number of points for curve integer npar_max parameter (npar_max = 100) real pi parameter (pi = 3.14159265) c arguments and common blocks integer npar, iflag double precision futil, chi2, par(*), grad(*) c local variables character*1 dummy character*80 infile / 'ingress_24_00.dat' / character*80 line character*80 title / ' ' / integer frame(max_sets) save frame integer i integer ierr integer ios integer istat integer j integer len integer lun integer set integer num_points(max_sets) save num_points integer num_sets save num_sets integer sym_type logical HEXIST logical store_plot / .false. / double precision c2(max_sets) save c2 double precision x(max_points, max_sets) double precision y(max_points,max_sets) double precision dx(max_points, max_sets) double precision dy(max_points, max_sets) double precision deriv(max_points, max_sets) save x, y, dx, dy, deriv real theta real sym_size / 0.1 / real xf (nf), yf (nf) real x_min, x_max real y_min, y_max real xs, ys, rs, xv, yv, rv real x4(max_points), y4(max_points) real dx4(max_points), dy4(max_points) save x4, y4, dx4, dy4 double precision xp, yp ! fit functions, defined below c begin if ( iflag .eq. 1 ) then ! read input data call KUPROS ('enter file name for input data', infile, len) lun = 30 open (unit=lun, file = infile, 1 form = 'formatted', status = 'old') set = 0 i = 0 ios = 0 do while ( ios .eq. 0 ) read (lun, fmt = '(a80)', iostat = ios) line if (ios .eq. 0 .and. .not. (line(1:1) .eq. '!') ) then if ( line(1:7) .eq. '# frame' ) then set = set + 1 i = 0 read(line(9:12), fmt='(i4)') frame(set) write (*, *) 'reading frame ', frame(set) else i = i + 1 read (line, *) x(i,set), y(i,set), & dx(i,set), dy(i,set), deriv(i,set) c write (*, *) 'reading values for i, set = ', i, set c write (*, *) x(i,set), y(i,set), c & dx(i,set), dy(i,set), deriv(i,set) num_points(set) = i endif endif end do close (lun) num_sets = set c check... c write (*, *) 'num_sets = ', num_sets c do j = 1, num_sets c write (*, *) 'frame ', frame(j), ' ...' c write (*, *) 'num_points = ', num_points(j) c write (*, *) 'x, y, dx, dy...' c do i = 1, num_points(j) c write (*, *) x(i,j), y(i,j), dx(i,j), dy(i,j) c end do c end do c define parameters c Define parameters alpha and beta, give initial values and step sizes. call MNPARM (1, 'rs', 335.d0, 0.1d0, 0.d0, 0.d0, ierr) call MNPARM (2, 'rv', 4.d0, 0.05d0, 0.d0, 0.d0, ierr) do i = 1, num_sets call MNPARM (i*4-1, 'xs', 335.d0, 1.d0, 0.d0, 0.d0, ierr) call MNPARM (i*4, 'ys', 201.d0, 1.d0, 0.d0, 0.d0, ierr) call MNPARM (i*4+1, 'xv', 48.d0, 0.5d0, 0.d0, 0.d0, ierr) call MNPARM (i*4+2, 'yv', 58.d0, 0.5d0, 0.d0, 0.d0, ierr) end do endif if ( iflag .eq. 2 ) then ! calculate gradient vector endif c calculate chi2 chi2 = 0. do j = 1, num_sets c2(j) = 0. do i = 1, num_points(j) if ( dx(i,j) .eq. 0. ) then c2(j) = c2(j) + & (y(i,j) - yp(x(i,j), y(i,j), par, j))**2 / dy(i,j)**2 elseif ( dy(i,j) .eq. 0. ) then c2(j) = c2(j) + & (x(i,j) - xp(x(i,j), y(i,j), par, j))**2 / dx(i,j)**2 endif end do chi2 = chi2 + c2(j) end do if ( iflag .eq. 3 ) then ! comes here to make plots call SET_PLOT_OPTIONS ! sets up sizes of various things call KUPROI ('store plot in ps file (0=no,1=yes)', store_plot) if ( store_plot ) call PS_OUTPUT c Set plot limits (hit carriage return to accept current values) x_min = 0. x_max = 120. y_min = 0. y_max = 120. title = 'pixels' c call KUPROR ('enter x_min', x_min) c call KUPROR ('enter x_max', x_max) c call KUPROR ('enter y_min', y_min) c call KUPROR ('enter y_max', y_max) c call KUPROS ('enter x axis title', title, len) c Set up the plot with HBOOK/HPLOT routines do set = 1, num_sets if ( HEXIST (1) ) call HDELET (1) title = 'frame ' write (title(7:11), fmt='(i5)') frame(set) call HBOOK1 (1, title, 2, x_min, x_max, 0.) call HMINIM (1, y_min) call HMAXIM (1, y_max) call HPLSET ('DMOD', 1.) ! solid line call HPLOT (1, ' ', 'HIST', 0) c plot the data values sym_type = 20 ! circle sym_size = 0.1 do i = 1, num_points(set) x4(i) = x(i,set) y4(i) = y(i,set) dx4(i) = dx(i,set) dy4(i) = dy(i,set) end do call HPLERR (x4, y4, dx4, dy4, num_points(set), ' ', & sym_type, sym_size) c plot the fitted function (evaluate at nf points) rs = par(1) rv = par(2) xs = par(4*set - 1) ys = par(4*set) xv = par(4*set + 1) yv = par(4*set + 2) do i = 1, nf theta = 2.*pi*FLOAT(i-1)/FLOAT(nf-1) xf(i) = xs + rs*COS(theta) yf(i) = ys + rs*SIN(theta) end do call HPLFUN (xf, yf, nf, ' ') do i = 1, nf theta = 2.*pi*FLOAT(i-1)/FLOAT(nf-1) xf(i) = xv + rv*COS(theta) yf(i) = yv + rv*SIN(theta) end do call HPLFUN (xf, yf, nf, ' ') c show parameter values and chi2 to left of plot call SHOW_PAR(c2, set, frame) call IGTERM ! force to screen write (*, *) 'Press any key to continue' read (*, fmt='(a1)') dummy end do ! 1, num_sets if ( store_plot ) call PS_OUTPUT ! closes ps file endif ! iflag .eq. 3 return END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function xp (x, y, par, set) c xp predicted x value of point on circle; take choice closest to x implicit NONE integer i, min_i, set double precision x, y double precision par (*) double precision xs, ys, rs, xv, yv, rv double precision d2, delta(4), pred(4), delta_min c begin rs = par(1) rv = par(2) xs = par(4*set-1) ys = par(4*set) xv = par(4*set+1) yv = par(4*set+2) d2 = rs**2 - (y - ys)**2 if ( d2 .gt. 0 ) then pred(1) = xs + sqrt(d2) pred(2) = xs - sqrt(d2) else pred(1) = 99999. pred(2) = 99999. endif d2 = rv**2 - (y - yv)**2 if ( d2 .gt. 0 ) then pred(3) = xv + sqrt(d2) pred(4) = xv - sqrt(d2) else pred(3) = 99999. pred(4) = 99999. endif do i = 1, 4 delta(i) = ABS(pred(i) - x) end do min_i = 1 delta_min = delta(1) do i = 2, 4 if ( delta(i) .lt. delta_min ) then min_i = i delta_min = delta(i) endif end do xp = pred(min_i) return END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function yp (x, y, par, set) c yp predicted y value of point on circle; take choice closest to y implicit NONE integer i, min_i, set double precision x, y double precision par (*) double precision xs, ys, rs, xv, yv, rv double precision d2, delta(4), pred(4), delta_min c begin rs = par(1) rv = par(2) xs = par(4*set-1) ys = par(4*set) xv = par(4*set+1) yv = par(4*set+2) d2 = rs**2 - (x - xs)**2 if ( d2 .gt. 0 ) then pred(1) = ys + sqrt(d2) pred(2) = ys - sqrt(d2) else pred(1) = 99999. pred(2) = 99999. endif d2 = rv**2 - (x - xv)**2 if ( d2 .gt. 0 ) then pred(3) = yv + sqrt(d2) pred(4) = yv - sqrt(d2) else pred(3) = 99999. pred(4) = 99999. endif do i = 1, 4 delta(i) = ABS(pred(i) - y) end do min_i = 1 delta_min = delta(1) do i = 2, 4 if ( delta(i) .lt. delta_min ) then min_i = i delta_min = delta(i) endif end do yp = pred(min_i) return END