program test_conic implicit none ! tests the conic section fitting routines ! declare character*1 data_type,conic_type character*40 string integer i,j,jend,n,nmax parameter (nmax=100) double precision x(nmax+3),y(nmax+3),x0,y0,angmin,angmax,ampl, & rmaj,rmin,tilt,xxx(nmax+3),yyy(nmax+3),sig(nmax+3) ! for the best fit geometry double precision xcenter,ycenter,tilt_out,semi_major,semi_minor,radius,flength, & eccen,xfoci1,yfoci1,xfoci2,yfoci2 ! for the results of the fitting the noisy data double precision a,asig,b,bsig,c,csig,d,dsig,e,esig,f,fsig,chi,coef(6) ! .popular formats 01 format(1x,a,1p10e12.4) 02 format(a,i3.3,a) 03 format(1x,a,1pe11.3,' +/-',1pe11.3) 04 format(1x,1p2e12.4) ! set the initial data type data_type = 'E' ! set the random number seed call random_seed(put=(/ 123456789, 987654321 /)) ! generate a noisy data set ampl = 0.03d0 x0 = 0.75d0 y0 = 0.25d0 tilt = 30.0d0 n = 40 if (data_type .eq. 'E') then rmaj = 1.25d0 rmin = 0.75d0 angmin = 0.0d0 angmax = 360.0d0 call generate_noisy_ellipse(x0,y0,rmaj,rmin,tilt,angmin,angmax,ampl,n,x,y,sig) else if (data_type .eq. 'H') then n = n - mod(n,2) rmaj = 0.30d0 rmin = 0.15d0 angmin = -2.0d0 angmax = 2.0d0 call generate_noisy_hyperbola(x0,y0,rmaj,rmin,tilt,angmin,angmax,ampl,n,x,y,sig) else if (data_type .eq. 'C') then rmaj = 0.75 angmin = 0.0d0 angmax = 360.0d0 call generate_noisy_circle(x0,y0,rmaj,angmin,angmax,ampl,n,x,y,sig) else if (data_type .eq. 'P') then rmaj = 0.2d0 angmin = 40.0d0 angmax = 90.0d0 call generate_noisy_parabola(x0,y0,rmaj,tilt,angmin,angmax,ampl,n,x,y,sig) else stop 'unknown conic section type' end if ! write out the noisy data write(string,02) 'noisy_data_',int(tilt),'.dat' open (unit=3,file=string,status='unknown') write(3,04) (x(i), y(i), i=1,n) close(unit=3) ! find the best fit call fitconic(x,y,n,a,asig,b,bsig,c,csig,d,dsig,e,esig,f,fsig,chi) ! report the fitting coefficients write(6,*) ' ' write(6,03) 'a=',a,asig,'b=',b,bsig,'c=',c,csig, & 'd=',d,dsig,'e=',e,esig,'f=',f,fsig write(6,01) 'chi squared:',chi ! now determine the conic section type coef = (/ a, b, c, d, e, f /) call check_type_of_conic(coef,conic_type) jend = min(361,nmax) jend = jend - mod(jend,2) ampl = 0.0d0 if (conic_type .eq. 'E') then call ellipse_geometry_from_coefs(coef,xcenter,ycenter,tilt_out, & semi_major,semi_minor,eccen,xfoci1,yfoci1,xfoci2,yfoci2) write(6,01) 'center :',xcenter,ycenter write(6,01) 'tilt :',tilt_out write(6,01) 'smajor axis :',semi_major write(6,01) 'sminor axis :',semi_minor write(6,01) 'eccentricity:',eccen write(6,01) 'foci1 at :',xfoci1,yfoci1 write(6,01) 'foci2 at :',xfoci2,yfoci2 xxx(1) = xcenter yyy(1) = ycenter xxx(2) = xfoci1 yyy(2) = yfoci1 xxx(3) = xfoci2 yyy(3) = yfoci2 angmin = 0.0d0 angmax = 360.0d0 call generate_noisy_ellipse(xcenter,ycenter,semi_major,semi_minor,tilt_out, & angmin,angmax,ampl,jend,xxx(4),yyy(4),sig(4)) else if (conic_type .eq. 'H') then call hyperbola_geometry_from_coefs(coef,xcenter,ycenter,tilt_out, & semi_major,semi_minor,eccen,xfoci1,yfoci1,xfoci2,yfoci2) write(6,01) 'center :',xcenter,ycenter write(6,01) 'tilt :',tilt_out write(6,01) 'smajor axis :',semi_major write(6,01) 'sminor axis :',semi_minor write(6,01) 'eccentricity:',eccen write(6,01) 'foci1 at :',xfoci1,yfoci1 write(6,01) 'foci2 at :',xfoci2,yfoci2 xxx(1) = xcenter yyy(1) = ycenter xxx(2) = xfoci1 yyy(2) = yfoci1 xxx(3) = xfoci2 yyy(3) = yfoci2 angmax = 2.0d0*acos(-1.0d0/eccen) angmin = -angmax write(6,*) jend call generate_noisy_hyperbola(xcenter,ycenter,semi_major,semi_minor,tilt_out, & angmin,angmax,ampl,jend,xxx(4),yyy(4),sig(4)) else if (conic_type .eq. 'C') then call ellipse_geometry_from_coefs(coef,xcenter,ycenter,tilt_out, & semi_major,semi_minor,eccen,xfoci1,yfoci1,xfoci2,yfoci2) write(6,01) 'center :',xcenter,ycenter write(6,01) 'radius :',semi_major xxx(1) = xcenter yyy(1) = ycenter xxx(2) = xcenter yyy(2) = ycenter xxx(3) = xcenter yyy(3) = ycenter angmin = 0.0d0 angmax = 360.0d0 call generate_noisy_circle(xcenter,ycenter,radius,angmin,angmax,ampl,jend,xxx(4),yyy(4),sig(4)) else if (conic_type .eq. 'P') then call parabola_geometry_from_coefs(coef,xcenter,ycenter,flength,tilt_out,xfoci1,yfoci1) write(6,01) 'center :',xcenter,ycenter write(6,01) 'tilt :',tilt_out write(6,01) 'flength :',flength write(6,01) 'foci1 at :',xfoci1,yfoci1 xxx(1) = xcenter yyy(1) = ycenter xxx(2) = xfoci1 yyy(2) = yfoci1 xxx(3) = xfoci1 yyy(3) = yfoci1 angmin = 20.0d0 angmax = 90.0d0 call generate_noisy_parabola(xcenter,ycenter,flength,tilt_out,angmin,angmax,ampl,jend,xxx(4),yyy(4),sig(4)) end if ! write out the best fit data write(string,02) 'curve_',int(tilt),'.dat' open (unit=3,file=string,status='unknown') write(3,04) (xxx(i), yyy(i), i=1,jend+3) close(unit=3) end subroutine fitconic(x,y,n, & a,asig,b,bsig,c,csig,d,dsig,e,esig,f,fsig,chisq) implicit none ! this routine fits data to a general conic ! input: ! x(1:n) = x-coordinates of the data ! y(1:n) = y coordinates of the data ! n = number of data points ! output: ! a,b,c,d,e,f = general conic equation coefficient ! asig to fsig = coefficient uncertainties, a +/- asig, b +/- bsig, and so on ! conic: a*x**2 + b*x*y + c*y**2 + d*x + e*y + 1 = 0 ! here we've normalized by the constant coefficient f, ! which will fail if the data pases through the origin (0,0) ! declare the pass integer n double precision x(n),y(n), & a,asig,b,bsig,c,csig,d,dsig,e,esig,f,fsig,chisq ! local variables integer i,nmax,npar parameter (nmax=500, npar=5) double precision u(nmax,npar),w(npar),v(npar,npar),bp(nmax), & ap(npar),cvm(npar,npar),wmax,thresh,sum,tol parameter (tol = 1.0d-12) ! bullet check if (n .gt. nmax) stop 'n > nmax in routine fitell' ! set the u matrix and right hand side ! 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 u(i,1) = x(i)**2 u(i,2) = x(i)*y(i) u(i,3) = y(i)**2 u(i,4) = x(i) u(i,5) = y(i) bp(i) = -1.0d0 enddo ! 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 ! return the uncertainties in the fitted parameters sqrt(cvm(i,i)) 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)**2 + ap(2)*y(i)**2 + ap(3)*x(i)*y(i) & + ap(4)*x(i) + ap(5)*y(i) chisq = chisq + (bp(i) - sum)**2 enddo chisq = chisq/(n - npar) ! set the return parameters a = ap(1) b = ap(2) c = ap(3) d = ap(4) e = ap(5) f = 1.0d0 ! and their uncertainties asig = sqrt(cvm(1,1)) bsig = sqrt(cvm(2,2)) csig = sqrt(cvm(3,3)) dsig = sqrt(cvm(4,4)) esig = sqrt(cvm(5,5)) fsig = 0.0d0 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 destroyed, 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 ! ! general utility routines below ! subroutine generate_noisy_ellipse(x0,y0,rmaj,rmin,tilt,angmin,angmax,ampl,n,x,y,sig) implicit none ! generates a noisy ellipse ! input: ! x0 = x coordinate of the ellipse center ! y0 = y coordinate of the ellipse center ! tilt = angle between the major axis and the x-axis in degrees ! rmaj = length of semi-major axis ! rmin = length of semi-minor axis ! angmin = minimum angle in degrees to start the noisy circle fro ! 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,rmaj,rmin,tilt,angmin,angmax,ampl,x(n),y(n),sig(n) ! local variables integer i double precision co,si,dangle,angle,xx,yy,rvec(2) ! go co = cosd(tilt) si = sind(tilt) dangle = (angmax - angmin)/float(n-1) do i = 1, n angle = angmin + dangle*dfloat(i-1) xx = rmaj * cosd(angle) yy = rmin * sind(angle) call gasdev(rvec) x(i) = co*xx - si*yy + x0 + ampl*rvec(1) y(i) = si*xx + co*yy + y0 + ampl*rvec(2) sig(i) = ampl if (ampl .eq. 0.0) sig(i) = 1.0d0 enddo return end subroutine generate_noisy_hyperbola(x0,y0,rmaj,rmin,tilt,angmin,angmax,ampl,n,x,y,sig) implicit none ! generates a noisy hyperbola ! input: ! x0 = x coordinate of the hyperbola center ! y0 = y coordinate of the hyperbola center ! tilt = angle between the major axis and the x-axis in degrees ! rmaj = length of semi-major axis ! rmin = length of semi-minor axis ! angmin = minimum angle in degrees to start the noisy circle fro ! 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 hyperbola ! n = number of data points to generate - should be an even number ! 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,rmaj,rmin,tilt,angmin,angmax,ampl,x(n),y(n),sig(n) ! local variables integer i,j double precision co,si,dangle,angle,xx,yy,rvec(2) ! put data on all four branches of the hyperbola co = cosd(tilt) si = sind(tilt) j = (n - mod(n,2))/2 dangle = (angmax - angmin)/float(j-1) do i = 1, j angle = angmin + dangle*dfloat(i-1) xx = rmaj * cosh(angle) yy = rmin * sinh(angle) call gasdev(rvec) x(i) = co*xx - si*yy + x0 + ampl*rvec(1) y(i) = si*xx + co*yy + y0 + ampl*rvec(2) sig(i) = ampl if (ampl .eq. 0.0) sig(i) = 1.0d0 xx = -xx yy = -yy call gasdev(rvec) x(i+j) = co*xx - si*yy + x0 + ampl*rvec(1) y(i+j) = si*xx + co*yy + y0 + ampl*rvec(2) sig(i+j) = ampl if (ampl .eq. 0.0) sig(i+j) = 1.0d0 enddo return end subroutine generate_noisy_parabola(x0,y0,flength,tilt,angmin,angmax,ampl,n,x,y,sig) implicit none ! generates a noisy parabola ! input: ! x0 = x coordinate of the parabola vertex ! y0 = y coordinate of the parabola vertex ! tilt = angle between the focal length and the x-axis in degrees ! flength = focal length ! angmin = minimum angle in degrees to start the noisy circle fro ! 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,flength,tilt,angmin,angmax,ampl,x(n),y(n),sig(n) ! local variables integer i,j double precision co,si,dangle,angle,rr2,xx,yy,rvec(2) ! put data on both branches of the parabola ! rr2 is the expression for a parabola in polar coordinates references from the origin ! the commented out rr2 is the parabolar in polars references from the foci co = cosd(tilt) si = sind(tilt) j = (n - mod(n,2))/2 dangle = (angmax - angmin)/float(j-1) do i=1,j angle = angmin + dangle*dfloat(i-1) rr2 = 4.0d0*flength*cosd(angle)/sind(angle)**2 ! rr2 = -2.0d0*flength/(1.0d0 + cosd(angle)) xx = rr2*cosd(angle) yy = rr2*sind(angle) call gasdev(rvec) x(j-i+1) = co*xx - si*yy + x0 + ampl*rvec(1) y(j-i+1) = si*xx + co*yy + y0 + ampl*rvec(2) sig(j-i+1) = ampl if (ampl .eq. 0.0) sig(i) = 1.0d0 ! different scatter on the other branch yy = -yy call gasdev(rvec) x(j+i) = co*xx - si*yy + x0 + ampl*rvec(1) y(j+i) = si*xx + co*yy + y0 + ampl*rvec(2) sig(j+i) = ampl if (ampl .eq. 0.0) sig(j+i) = 1.0d0 enddo return 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 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. fac = sqrt(-2.0d0*log(rsq)/rsq) rvec(1) = v1*fac rvec(2) = v2*fac return end subroutine ellipse_geometry_from_coefs(coef, & x0,y0,tilt,rmaj,rmin,eccen,xfoci1,yfoci1,xfoci2,yfoci2) 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 ellipse center ! y0 = y coordinate of the ellipse center ! tilt = angle between the major axis and the x-axis ! rmaj = length of semi-major axis ! rmin = length of semi-minor axis ! eccen = eccentricity ! xfoci1 = x coordinate of foci1 ! yfoci1 = y coordinate of foci1 ! xfoci2 = x coordinate of foci2 ! yfoci2 = y coordinate of foci2 ! declare the pass double precision coef(6),x0,y0,tilt,rmaj,rmin,eccen,xfoci1,yfoci1,xfoci2,yfoci2 ! local variables double precision a,b,c,d,e,f,co,si,ar,cr,dr,er,fr,bigj,fff,t1,t2,foci_dist ! popular format statements 10 format(1x,1p10e14.6) ! put the coefficients into their classical form a = coef(1) b = coef(2) c = coef(3) d = coef(4) e = coef(5) f = coef(6) ! tilt - we may be off by 90 degrees but we'll adjust later tilt = 0.5d0 * atand(b/(a - c)) co = cosd(tilt) si = sind(tilt) ! conic coefficients in the rotated frame ar = a*co*co + b*si*co + c*si*si cr = a*si*si - b*si*co + c*co*co dr = d*co + e*si er = -d*si + e*co fr = f ! center ! this can fail if the data is exactly a parabola bigj = b*b - 4.0d0*a*c x0 = (2.0d0*c*d - b*e) / bigj y0 = (2.0d0*a*e - b*d) / bigj ! semi-major and semi-minor axes fff = (c*d*d + a*e*e - b*d*e + b*b*f - 4.0d0*a*c*f)/bigj t1 = sqrt(abs(fff/ar)) t2 = sqrt(abs(fff/cr)) ! if the tilt angle was correct if (t1 .gt. t2) then rmaj = t1 rmin = t2 ! if the tilt angle is not correct, swap the major and minor axes ! and adjust the tilt else rmaj = t2 rmin = t1 tilt = tilt + 90.0d0 co = cosd(tilt) si = sind(tilt) end if ! distance from center to the foci and the eccentricity foci_dist = sqrt(rmaj**2 - rmin**2) eccen = foci_dist/rmaj ! coordinates of the two foci xfoci1 = x0 + co*foci_dist yfoci1 = y0 + si*foci_dist xfoci2 = x0 - co*foci_dist yfoci2 = y0 - si*foci_dist return end subroutine hyperbola_geometry_from_coefs(coef, & x0,y0,tilt,rmaj,rmin,eccen,xfoci1,yfoci1,xfoci2,yfoci2) 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 hyperbola center ! y0 = y coordinate of the hyperbola center ! tilt = angle between the major axis and the x-axis ! rmaj = length of semi-major axis ! rmin = length of semi-minor axis ! eccen = eccentricity ! xfoci1 = x coordinate of foci1 ! yfoci1 = y coordinate of foci1 ! xfoci2 = x coordinate of foci2 ! yfoci2 = y coordinate of foci2 ! declare the pass double precision coef(6),x0,y0,tilt,rmaj,rmin,eccen,xfoci1,yfoci1,xfoci2,yfoci2 ! local variables double precision a,b,c,d,e,f,co,si,ar,cr,dr,er,fr,bigj,fff,t1,t2,foci_dist ! put the coefficients into their classical form a = coef(1) b = coef(2) c = coef(3) d = coef(4) e = coef(5) f = coef(6) ! tilt - we may be off by 90 degrees but we'll adjust later tilt = 0.5d0 * atand(b/(a - c)) co = cosd(tilt) si = sind(tilt) ! conic coefficients in the rotated frame ar = a*co*co + b*si*co + c*si*si cr = a*si*si - b*si*co + c*co*co dr = d*co + e*si er = -d*si + e*co fr = f ! center ! this can fail if the data is exactly a parabola bigj = b*b - 4.0d0*a*c x0 = (2.0d0*c*d - b*e) / bigj y0 = (2.0d0*a*e - b*d) / bigj ! semi-major and semi-minor axes fff = (c*d*d + a*e*e - b*d*e + b*b*f - 4.0d0*a*c*f)/bigj t1 = sqrt(abs(fff/ar)) t2 = sqrt(abs(fff/cr)) ! if the tilt angle was correct if (t1 .gt. t2) then rmaj = t1 rmin = t2 ! if the tilt angle is not correct, swap the major and minor axes ! and adjust the tilt else rmaj = t2 rmin = t1 tilt = tilt + 90.0d0 co = cosd(tilt) si = sind(tilt) end if ! distance from center to the foci and the eccentricity foci_dist = sqrt(rmaj**2 + rmin**2) eccen = foci_dist/rmaj ! coordinates of the two foci xfoci1 = x0 + co*foci_dist yfoci1 = y0 + si*foci_dist xfoci2 = x0 - co*foci_dist yfoci2 = y0 - si*foci_dist 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 parabola_geometry_from_coefs(coef,x0,y0,flength,tilt,xfoci1,yfoci1) 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 ! flength = focal length ! tilt = angle between x-axis and axis of of symmetry ! xfoci1 = x coordinate of the foci ! yfoci1 = y coordinate of the foci ! declare the pass double precision coef(6),x0,y0,flength,tilt,xfoci1,yfoci1 ! local variables double precision a,b,c,d,e,f,co,si,bprime,dprime,eprime,t1,t2,large parameter (large=1.0d6) ! put the coefficients into their classical form a = coef(1) b = coef(2) c = coef(3) d = coef(4) e = coef(5) f = coef(6) !..trial value for the tilt, it may or may not be off by +/- pi tilt = 0.5d0 * atand(b/(a - c)) co = cosd(tilt) si = sind(tilt) bprime = a*si**2 - b*co*si + c*co**2 dprime = d*co + e*si flength = -0.25d0 * dprime/bprime !..if the focal length is monstrous, we are off by pi/2 100 if (flength .gt. large) then tilt = tilt + 90.0d0 co = cosd(tilt) si = sind(tilt) bprime = a*si**2 + c*co**2 - b*co*si dprime = d*co + e*si flength = -0.25d0 * dprime/bprime end if !..if the focal length is negative, we are off by pi if (flength .lt. 0.0) then tilt = tilt + 180.0d0 co = cosd(tilt) si = sind(tilt) bprime = a*si**2 + c*co**2 - b*co*si dprime = d*co + e*si flength = -0.25d0 * dprime/bprime goto 100 end if if (tilt .gt. 360.0) tilt = tilt - 360.0 !..carry on eprime = e*co - d*si t2 = -0.5d0 * eprime/bprime t1 = -(1.0d0 + eprime*t2 + bprime*t2**2)/dprime x0 = t1*co - t2*si y0 = t1*si + t2*co xfoci1 = x0 + flength*co yfoci1 = y0 + flength*si return end subroutine check_type_of_conic(coef,conic_type) 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 character*1 conic_type double precision coef(*) ! local variables double precision a,b,c,d,e,f,m11,m23,det,small parameter (small = 1.0d-10) ! initialize and map the vector conic_type = ' ' 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 ((abs(a-c) .le. small) .and. (abs(b) .le. small)) then write(6,*) 'circle' conic_type = 'C' else write(6,*) 'ellipse' conic_type = 'E' end if endif else if (m11 .lt. -small) then write(6,*) 'hyperbola' conic_type = 'H' else write(6,*) 'parabola' conic_type = 'P' end if endif return end