program fit_circle implicit none ! tests the circle specific fitting ! for generating and storing the noisy data integer i,n,nmax parameter (nmax=100) double precision x(nmax+1),y(nmax+1),sig(nmax+1),x0,y0,angmin,angmax, & dangle,angle,xx,yy,ampl,radius,rr,rvec(2) ! for the circle fitter integer iprint double precision xcenter,ycenter,rradius,coef(6),ucoef(6),chisq ! for plotting the fitted conic character*40 string integer jend double precision xxx(nmax+3),yyy(nmax+3) ! popular formats 01 format(1x,a,1p10e12.4) 03 format(1x,1pe11.3,' +/- ',1pe11.3,' * x**2 +'/, & 1x,1pe11.3,' +/- ',1pe11.3,' * y**2 +',/, & 1x,1pe11.3,' +/- ',1pe11.3,' * x*y +',/, & 1x,1pe11.3,' +/- ',1pe11.3,' * x +',/, & 1x,1pe11.3,' +/- ',1pe11.3,' * y +',/, & 1x,1pe11.3,' +/- ',1pe11.3,' = 0') 04 format(1x,1p2e12.4) ! get some random numbers call random_seed(put=(/ 123456789, 987654321 /)) ! generate a noisy data set x0 = 0.75d0 y0 = 0.25d0 radius = 1.25d0 angmin = 0.0d0 angmax = 360.0d0 ampl = 0.03d0 n = 40 call generate_noisy_circle(x0,y0,radius,angmin,angmax,ampl,n,x,y,sig) ! fit a cricle call fitcrc(x, y, sig, n, & xcenter,ycenter,rradius,coef,ucoef,chisq) ! report the results write(6,*) 'geometric results of fitting the data to an ellipse:' write(6,01) 'center at :',xcenter,ycenter write(6,01) 'radius :',rradius write(6,*) ' ' write(6,01) 'equation of best fit conic:' write(6,03) (coef(i),ucoef(i), i=1,6) ! store the fitted curve, put the center as the first point xxx(1) = xcenter yyy(1) = ycenter jend = min(361,nmax) angmin = 0.0d0 angmax = 360.0d0 ampl = 0.0d0 call generate_noisy_circle(xcenter,ycenter,rradius,angmin,angmax,ampl,jend,xxx(2),yyy(2),sig) ! write out the random points and the best fit for plotting write(string,67) 'noisy_data_circle01.dat' 67 format(a) open (unit=3,file=string,status='unknown') write(3,04) (x(i), y(i), i=1,n) close(unit=3) write(string,67) 'curve_circle01.dat' open (unit=3,file=string,status='unknown') write(3,04) (xxx(i), yyy(i), i=1,jend+1) close(unit=3) end subroutine generate_noisy_circle(x0,y0,radius,angmin,angmax,ampl,n,x,y,sig) implicit none ! generates a noisy circle ! input: ! x0 = x coordinate of the center ! y0 = y coordinate of the center ! r = radius ! angmin = minimum angle in degrees to start the noisy circle from ! angmax = maximim angle in degrees to end the noisy circle from ! ampl = amplitude of the gaussian noie to be added - set to 0.0 for a perfect ellipse ! n = number of data points to generate ! output ! x = x-coordinate of data ! y = y-coordinate of data ! sig = gaussian uncertainty in the data ! declare the pass integer n double precision x0,y0,radius,angmin,angmax,ampl,x(n),y(n),sig(n) ! local variables integer i double precision dangle,angle,rr,rvec(2) ! go dangle = (angmax - angmin)/float(n-1) do i = 1, n angle = angmin + dangle*dfloat(i-1) call gasdev(rvec) x(i) = radius*cosd(angle) + x0 + ampl*rvec(1) y(i) = radius*sind(angle) + y0 + ampl*rvec(2) sig(i) = ampl if (ampl .eq. 0.0) sig(i) = 1.0d0 enddo return end subroutine fitcrc(x,y,sig,n,x0,y0,radius,coef,ucoef,chisq) implicit none ! this routine fits data to a circle. ! input: ! x(1:n) = x-coordinates of the data ! y(1:n) = y coordinates of the data ! sig(i) = uncertainty in the coordinate positions ! n = number of data points ! output: ! x0 = x coordinate of the circle center ! y0 = y coordinate of the circle center ! radius = radius of the circle ! coef = coefficients of best fit conic ! ucoef = uncertainty in coefficients of best fit conic ! chisq = chi-square per degree of freedom ! circle equation: x**2 + y**2 + d*x + e*y + f = 0 ! declare the pass integer n double precision x(n),y(n),sig(n),x0,y0,radius,coef(6),ucoef(6),chisq ! local variables integer i,j,nmax,npar parameter (nmax=500, npar=3) double precision u(nmax,npar),w(npar),v(npar,npar),bp(nmax), & ap(npar),cvm(npar,npar),wmax,thresh,sum,r2,tmp,tol parameter (tol = 1.0d-6) ! bullet check if (n .gt. nmax) stop 'n > nmax in routine fitcrc' ! set the u matrix and right hand side bp ! the entries of u are simply the terms in front of the parameters to be fit ! and the right hand side is everything else moved across the equal sign do i = 1,n tmp = 1.0d0/sig(i) u(i,1) = x(i) * tmp u(i,2) = y(i) * tmp u(i,3) = 1.0d0 * tmp bp(i) = -(x(i)**2 + y(i)**2) * tmp enddo ! do a singular value decomposition call svdcmp(u,n,npar,nmax,npar,w,v) !massage the singular values wmax = 0.0d0 do i=1,npar if (w(i) .gt. wmax) wmax = w(i) enddo thresh = tol * wmax do i=1,npar if (w(i) .lt. thresh) w(i) = 0.0d0 enddo ! backsubstitution and covariance matrix ! calls to double precision versions of the numerical recipes routines call svbksb(u,w,v,n,npar,nmax,npar,bp,ap) call svdvar(v,npar,npar,w,cvm,npar) ! evaluate the chi-square per degree of freedom chisq = 0.0d0 do i=1,n sum = ap(1)*x(i) + ap(2)*y(i) + ap(3) chisq = chisq + ((bp(i) - sum)/sig(i))**2 enddo chisq = chisq/(n - npar) ! set the return parameters ! the circle equation coefficients coef(1) = 1.0d0 coef(2) = 0.0d0 coef(3) = 1.0d0 coef(4) = ap(1) coef(5) = ap(2) coef(6) = ap(3) ! and their standard deviation uncertainties ucoef(1) = 0.0d0 ucoef(2) = 0.0d0 ucoef(3) = 0.0d0 ucoef(4) = sqrt(cvm(1,1)) ucoef(5) = sqrt(cvm(2,2)) ucoef(6) = sqrt(cvm(3,3)) ! get the circle geometry call circle_geometry_from_coefs(coef,x0,y0,radius) return end subroutine circle_geometry_from_coefs(coef,x0,y0,radius) implicit none ! input: ! coef = the six coeficients to the conic equation ! coef(1)*x**2 + coef(2)*y**2 + coef(3)*x*y + coef(4)*x + coef(5)*y + coef(6) = 0 ! ! output: ! x0 = x coordinate of the circle center ! y0 = y coordinate of the circle center ! radius = radius of the circle ! declare the pass double precision coef(6),x0,y0,radius ! local variables double precision r2 ! circle center and radius x0 = -0.5d0 * coef(4) y0 = -0.5d0 * coef(5) r2 = x0*x0 + y0*y0 - coef(6) radius = 1.0d30 if (r2 .gt. 0.0d0) radius = sqrt(r2) return end subroutine svdcmp(a,m,n,mp,np,w,v) implicit none ! given a matrix a, logical dimensions m by n, physical dimension mp by np, ! this routine computes its singular value decomposition ! a = u * w * vtranspose ! the matrix u replaces a on output. The diagonal of singular values ! is output as the vector w. the matrix v, not its transpose, is also output. ! m must be greater than or equal to n; if it is smaller, then a should ! be filled up to a square with zero rows ! declare integer m,n,mp,np,i,j,k,l,its,nmax,nm,jj parameter (nmax = 500) double precision a(mp,np),w(np),v(np,np),rv1(nmax),scale, & anorm,c,f,g,h,s,x,y,z,pythag ! householder reduction to bidiagonal form g = 0.0d0 scale = 0.0d0 anorm = 0.0d0 do i=1,n l = i + 1 rv1(i) = scale*g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if (i .le. m) then do k=i,m scale = scale + abs(a(k,i)) enddo if (scale .ne. 0.0) then do k=i,m a(k,i) = a(k,i)/scale s = s + a(k,i)*a(k,i) enddo f = a(i,i) g = -sign(sqrt(s),f) h = f*g - s a(i,i) = f - g do j=l,n s = 0.0d0 do k=i,m s = s + a(k,i)*a(k,j) enddo f = s/h do k=i,m a(k,j) = a(k,j) + f*a(k,i) enddo enddo do k=i,m a(k,i) = scale*a(k,i) enddo end if end if w(i) = scale *g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if ((i .le. m) .and. (i .ne. n)) then do k=l,n scale = scale + abs(a(i,k)) enddo if (scale .ne. 0.0) then do k=l,n a(i,k) = a(i,k)/scale s = s + a(i,k)*a(i,k) enddo f = a(i,l) g = -sign(sqrt(s),f) h = f*g - s a(i,l) = f - g do k=l,n rv1(k) = a(i,k)/h enddo do j=l,m s = 0.0d0 do k=l,n s = s + a(j,k)*a(i,k) enddo do k=l,n a(j,k) = a(j,k) + s*rv1(k) enddo enddo do k=l,n a(i,k) = scale*a(i,k) enddo end if end if anorm = max(anorm,(abs(w(i))+abs(rv1(i)))) ! end of householder reduction to bidiagonal form enddo ! accumulation of the right hand transformations do i=n,1,-1 if (i .lt. n) then ! double division to avoid underflow if (g .ne. 0.0) then do j=l,n v(j,i) = (a(i,j)/a(i,l))/g enddo do j=l,n s = 0.0d0 do k=l,n s = s + a(i,k)*v(k,j) enddo do k=l,n v(k,j) = v(k,j) + s*v(k,i) enddo enddo end if do j=l,n v(i,j) = 0.0d0 v(j,i) = 0.0d0 enddo end if v(i,i) = 1.0d0 g = rv1(i) l = i ! end accumulation of the right hand transformations enddo ! accumulate the left hand transforms do i=min(m,n),1,-1 l = i+1 g = w(i) do j=l,n a(i,j) = 0.0d0 enddo if (g .ne. 0.0) then g=1.0d0/g do j=l,n s = 0.0d0 do k=l,m s = s + a(k,i)*a(k,j) enddo f = (s/a(i,i))*g do k=i,m a(k,j) = a(k,j) + f*a(k,i) enddo enddo do j=i,m a(j,i) = a(j,i)*g enddo else do j= i,m a(j,i) = 0.0d0 enddo end if a(i,i) = a(i,i) + 1.0d0 ! end accumulate the left hand transforms enddo ! diagonalization of the bidiagonal form ! loop over singular values do k=n,1,-1 ! loop over allowed iterations do its=1,30 ! test for splitting is here, note rv(1) is always zero do l=k,1,-1 nm = l-1 if ((abs(rv1(l))+anorm) .eq. anorm) goto 2 if ((abs(w(nm))+anorm) .eq. anorm) goto 1 enddo ! cacellation of rv1(l) if l .gt. 1 1 c = 0.0d0 s = 1.0d0 do i=l,k f = s * rv1(i) rv1(i) = c * rv1(i) if ((abs(f)+anorm) .eq. anorm) goto 2 g = w(i) h = pythag(f,g) w(i) = h h = 1.0d0/h c = (g*h) s = -(f*h) do j=1,m y = a(j,nm) z = a(j,i) a(j,nm) = (y*c) + (z*s) a(j,i) = -(y*s) + (z*c) enddo enddo ! convergence check 2 z = w(k) if (l .eq. k) then if (z .lt. 0.0) then w(k) = -z do j=1,n v(j,k) = -v(j,k) enddo end if goto 3 end if if (its .eq. 30) stop 'no convergence in 30 iters in svdcmp' ! shift the bottom 2 by 2 minor x = w(l) nm = k-1 y = w(nm) g = rv1(nm) h = rv1(k) f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0d0*h*y) g = pythag(f,1.0d0) f = ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x ! now the qr transformation c = 1.0d0 s = 1.0d0 do j=l,nm i = j+1 g = rv1(i) y = w(i) h = s*g g = c*g z = pythag(f,h) rv1(j) = z c = f/z s = h/z f = (x*c) + (g*s) g = -(x*s) + (g*c) h = y*s y = y*c do jj=1,n x = v(jj,j) z = v(jj,i) v(jj,j) = (x*c)+(z*s) v(jj,i) = -(x*s)+(z*c) enddo z = pythag(f,h) w(j) = z ! the rotation can be arbitrary if z is zero if (z .ne. 0.0) then z = 1.0d0/z c = f*z s = h*z end if f = (c*g)+(s*y) x = -(s*g)+(c*y) do jj=1,m y = a(jj,j) z = a(jj,i) a(jj,j) = (y*c)+(z*s) a(jj,i) =-(y*s)+(z*c) enddo ! end of qr transform enddo rv1(l) = 0.0d0 rv1(k) = f w(k) = x ! end of iteration loop enddo 3 continue ! end of loop over singular values enddo return end double precision function pythag(a,b) implicit none ! computes sqrt(a**2 + b**2) without destructive underflow or overflow ! declare double precision a,b,absa,absb absa = abs(a) absb = abs(b) if (absa .gt. absb) then pythag = absa * sqrt(1.0d0 + (absb/absa)**2) else if (absb .eq. 0) then pythag = 0.0d0 else pythag = absb * sqrt(1.0d0 + (absa/absb)**2) end if end if return end subroutine svbksb(u,w,v,m,n,mp,np,b,x) implicit none ! solves ax=b for x when a is specified by the arrays u, w and v as returned ! by routine svdcmp. a is of logical dimension m by n and physical dimension ! mp by np. b is the input right hand side. x on output is the solution ! vector. no input quantities are destoyed, so one may call sequentially ! with different b's ! declare integer m,n,mp,np,nmax,i,j,jj parameter (nmax = 500) double precision u(mp,np),w(np),v(np,np),b(mp),x(np),tmp(nmax),s ! compute utranspose * b do j=1,n s = 0.0d0 ! nonzero result only if wj is nonzero; and divide if (w(j) .ne. 0.0) then do i=1,m s = s + u(i,j)*b(i) enddo s = s/w(j) end if tmp(j) = s enddo ! and matrix multiply to get the answer do j=1,n s = 0.0d0 do jj=1,n s = s + v(j,jj)*tmp(jj) enddo x(j) = s enddo return end subroutine svdvar(v,ma,np,w,cvm,ncvm) implicit none ! evaluates the covariance matrix cvm of the fit for ma parameters obtained ! by the routine svdfit, call this routine with matrices v and w as returned ! by svdfit. np and ncvm are the physical dimensions of v,w and cvm. ! declare integer ma,np,ncvm,mmax,i,j,k parameter (mmax=20) double precision v(np,np),w(np),cvm(ncvm,ncvm),wti(mmax),sum ! initialize do i=1,ma wti(i) = 0.0d0 if (w(i) .ne. 0.0) wti(i) = 1.0d0/(w(i)*w(i)) enddo ! sum the contributions to the covariant matrix; recipe equation 14.3.20 do i=1,ma do j=1,i sum = 0.0d0 do k=1,ma sum = sum + v(i,k)*v(j,k)*wti(k) enddo cvm(i,j) = sum cvm(j,i) = sum enddo enddo return end subroutine check_type_of_conic(coef) implicit none ! figure out which conic section it appears to be ! a*x**2 + b*y**2 + c*x*y + d*x + e*y + f = 0 ! http://mathworld.wolfram.com/QuadraticCurve.html ! declare the pass double precision coef(*) ! local variables double precision a,b,c,d,e,f,m11,m23,det,small parameter (small = 1.0d-10) ! map the vactor a = coef(1) b = coef(2) c = coef(3) d = coef(4) e = coef(5) f = coef(6) ! | F D/2 E/2 | ! M = | D/2 A B/2 | ! | E/2 B/2 C | m11 = a*c - 0.25d0*b*b m23 = 0.25d0*d*e - 0.5d0*b*f det = a*c*f + 0.25d0*(b*d*e - b*b*f - c*d*d - a*e*e) ! conic is degenerate if (det .eq. 0.0) then if (m11 .gt. small) then write(6,*) 'point' else if (m11 .lt. -small) then if (a + c .eq. 0.0) then write(6,*) 'perpendicular lines' else write(6,*) 'intersecting lines' end if else if (m23 .eq. 0.0d0) then write(6,*) 'single line' else write(6,*) 'two parallel lines' end if end if ! conic is an ellipse, hyperbola or parabola else if (m11 .gt. small) then if (a*det .gt. 0.0d0) then write(6,*) 'nothing' else if ((a .eq. c) .and. (b .eq. 0.0d0)) then write(6,*) 'circle' else write(6,*) 'ellipse' end if endif else if (m11 .lt. -small) then write(6,*) 'hyperbola' else write(6,*) 'parabola' end if endif return end subroutine gasdev(rvec) implicit none ! returns two gaussian distributed deviates with zero mean and unit variance ! declare double precision harvest,v1,v2,fac,gset,rsq,rvec(2) ! here we go 1 call random_number(rvec) v1 = 2.0d0 * rvec(1) - 1.0d0 v2 = 2.0d0 * rvec(2) - 1.0d0 ! see if they are in the unit circle, if not try again rsq = v1**2 + v2**2 if (rsq .ge. 1.0 .or. rsq .eq. 0.0) goto 1 ! now make the gaussian transform, recipe equation 7.2.10, to get two ! random deviates. return one and save the other for the next time ! then set the flag fac = sqrt(-2.0d0*log(rsq)/rsq) rvec(1) = v1*fac rvec(2) = v2*fac return end