C Test finding size of gaps between 3 interior circles parameter ( mmax=105,M2=MMAX+MMAX ) real x(mmax), y(mmax), r(mmax), XX(M2) CHARACTER FNAM*11 FNAM = 'UCIC00.TXT' C Name and version of program write (*,*) ' Program CIREXC, Version Date: December 5 2005' C M = 0 10 continue WRITE (*,*) 'N:' READ (*,*,end=10,err=10) N IF ( N .gt. 99 ) goto 10 IF ( N .lt. 0 ) then M = -N write (*,*) ' Minimum Circle:', M+1, ' N:' read (*,*) N if ( N .le. M ) goto 10 FNAM = 'B_00_00.txt' WRITE (FNAM(3:4),1002) N WRITE (FNAM(6:7),1002) M else WRITE (FNAM(5:6),1002) N endif 1002 FORMAT ( I2.2) OPEN ( UNIT=10, FORM='FORMATTED', FILE=FNAM, & STATUS='OLD', IOSTAT=IOS ) IF ( IOS .NE. 0 ) STOP 'Error opening file' REWIND ( UNIT=10 ) READ (10,*) K, DBEST IF ( N .NE. K ) STOP 'NON-FITTING CIRCLE NUMBER' C READ COORDINATES FROM FILE (NORMALIZED POSITION IS ASSUMED) L = 0 if ( m .eq. 0 ) then DO 5 K = 1, N READ (10,*) I, X(K), Y(K) R(K) = DBEST * FLOAT ( K ) L = L + 1 XX(L) = X(K) L = L + 1 XX(L) = Y(K) 5 CONTINUE else DBEST = 1.0 / DBEST DO 6 K = M+1, N READ (10,*) I, X(K), Y(K) X(K) = X(K) * DBEST Y(K) = Y(K) * DBEST R(K) = FLOAT ( K ) * DBEST L = L + 1 XX(L) = X(K) L = L + 1 XX(L) = Y(K) 6 CONTINUE endif C close ( unit=10 ) C 100 continue C Loops over all combinations of 3 circles do 110 I1 = M+1, N-2 do 120 I2 = I1+1, N-1 do 130 I3 = I2+1, N call c3cent ( N, X, Y, R, I1, I2, I3, XC, YC, RC ) C write (*,*) ' Center: ', XC, YC C write (*,*) ' Radius: ', RC, '. scaled: ', RC/DBEST C Check if circle remains inside container if ( xc**2 + yc**2 + RC**2 .gt. 1.0 ) then C Diagnostic output C write (*,*) ' Elmination outside:', C & sqrt ( xc**2 + yc**2 + RC**2 ) goto 130 endif C Check if other circles are cut or fall into the interior of the C gap-filling circle do 200 i = M+1, n if ( i .eq. i1 .or. i .eq. i2 .or. i .eq. i3 ) goto 200 if ( sqrt ((x(i)-xc)**2 + (y(i)-yc)**2) & .lt. R(i)+RC ) goto 130 200 continue write (*,2001) i1, i2, i3, xc, yc, rc/DBEST 2001 format ( 3I3, 2F12.6, F12.4 ) 130 continue 120 continue 110 continue 300 continue write (*,*) ' 3 circle indices' read (*,*) i1, i2, i3 if ( i1 .eq. 0 ) goto 400 call c3cent ( N, X, Y, R, I1, I2, I3, XC, YC, RC ) write (*,*) ' Center: ', XC, YC write (*,*) ' Radius: ', RC, '. scaled: ', RC/DBEST if ( xc**2 + yc**2 + RC**2 .gt. 1.0 ) then C Diagnostic output write (*,*) ' Elmination outside:', & sqrt ( xc**2 + yc**2 + RC**2 ) endif do 310 i = M+1, n if ( i .eq. i1 .or. i .eq. i2 .or. i .eq. i3 ) goto 310 if ( sqrt ((x(i)-xc)**2 + (y(i)-yc)**2) & .lt. R(i)+RC ) then write (*,*) ' Collision' write (*,*) i, sqrt ((x(i)-xc)**2 + (y(i)-yc)**2), R(I)+RC endif 310 continue goto 300 400 continue C write (*,*) ' 2 circle indices' C read (*,*) i1, i2 C if ( i1 .eq. 0 ) goto 999 C Loops over all combinations of 2 circles do 500 i1 = M+1, N-1 do 510 i2 = i1+1, N call c20cen ( N, X, Y, R, I1, I2, XC1, YC1, RC1, & XC2, YC2, RC2 ) C Paranoia: Check for crossing of container circle if ( xc1**2 + yc1**2 + RC1**2 .gt. 1.0 ) then C Diagnostic output write (*,*) ' Elmination outside:', & sqrt ( xc1**2 + yc1**2 + RC1**2 ) endif if ( xc2**2 + yc2**2 + RC2**2 .gt. 1.0 ) then C Diagnostic output write (*,*) ' Elmination outside:', & sqrt ( xc2**2 + yc2**2 + RC2**2 ) endif C C Collisions with other circles C do 409 irep = 1, 2 C Count collisions during first pass; if more than 1 collision is C found, skip second pass, in which printing is done if ( irep .eq. 1 ) then ncolis = 0 else if ( ncolis .gt. 2 ) goto 411 C write (*,*) ' Center1: ', XC1, YC1 C write (*,*) ' Radius1: ', RC1, '. scaled: ', RC1/DBEST write (*,2001) i1, i2, -1, xc1, yc1, rc1/DBEST endif do 410 i = M+1, n if ( i .eq. i1 .or. i .eq. i2 ) goto 410 if ( sqrt ((x(i)-xc1)**2 + (y(i)-yc1)**2) & .lt. R(i)+RC1 ) then ncolis = ncolis + 1 if ( irep .eq. 2 ) then write (*,*) ' Collision: ', & i, sqrt ((x(i)-xc1)**2 + (y(i)-yc1)**2), R(I)+RC1 endif endif 410 continue 409 continue 411 continue C do 419 irep = 1, 2 if ( irep .eq. 1 ) then ncolis = 0 else if ( ncolis .gt. 1 ) goto 421 C write (*,*) ' Center2: ', XC2, YC2 C write (*,*) ' Radius2: ', RC2, '. scaled: ', RC2/DBEST write (*,2001) i1, i2, -2, xc2, yc2, rc2/DBEST endif do 420 i = M+1, n if ( i .eq. i1 .or. i .eq. i2 ) goto 420 if ( sqrt ((x(i)-xc2)**2 + (y(i)-yc2)**2) & .lt. R(i)+RC2 ) then ncolis = ncolis + 1 if ( irep .eq. 2 ) then write (*,*) ' Collision: ', & i, sqrt ((x(i)-xc2)**2 + (y(i)-yc2)**2), R(I)+RC2 endif endif 420 continue 419 continue 421 continue C goto 400 510 continue 500 continue 999 continue if ( M .gt. 0 ) then C write preliminary rescaled result file FNAM = 'WCIC00.TXT' WRITE (FNAM(5:6),1002) N OPEN ( UNIT=10, FORM='FORMATTED', FILE=FNAM, & STATUS='UNKNOWN', IOSTAT=IOS ) IF ( IOS .NE. 0 ) STOP 'Error opening file' write (10,*) n, dbest, 1.0/dbest do 600 i = M+1, N write (10,*) I, x(i), y(i) 600 continue endif end C------------------------------------------------------------- subroutine c3cent ( M, X, Y, R, I1, I2, I3, XC, YC, RC ) C Determine center and radius of the circle touching 3 circles, C such all 3 circles are outside the central circle C In the case that the 3 other circles are mutually tangent, C this is known as the inner Soddy circle C real x(M), Y(M), R(M) real P(2), EPS(2) real x1, x2, x3, y1, y2, y3, r1, r2, r3 common /coor3c/ x1, x2, x3, y1, y2, y3, r1, r2, r3 external dist3 x1 = x(i1) x2 = x(i2) x3 = x(i3) y1 = y(i1) y2 = y(i2) y3 = y(i3) r1 = r(i1) r2 = r(i2) r3 = r(i3) C C xc = (((y1 - y3)*(y2 - y3) + x3**2)*(y1 - y2) - (y1 - y3) C & *x2**2 + (y2 - y3)*x1**2 + r3**2*y2 + (y1 - y3)*r2**2 C & - (y2 - y3)*r1**2 - r3**y1) / C & (2*((y1 - y2)*x3 - (y1 - y3)*x2 + (y2 - y3)*x1)) C C yc = (((y2 + y3)*(y2 - y3) - x3**2 + x2**2 - (x2 - x3)* C & x1)*x1 + ((y1 + y2)*(y1 - y2) - x2**2)*x3 - ((y1 + C & y3)*(y1 - y3) - x3**2)*x2 + (x1 - x2)*r3**2 - C & (x1 - x3)*r2**2 + (x2 - x3)*r1**2) / C & (2*((y1 - y2)* x3 - (y1 - y3)*x2 + (y2 - y3)*x1)) C C Until a correct formula becomes available, solve minimization C problem "Make squared radius differences as small as possible" EPS(1) = 0.0 EPS(2)= 0.0 d12 = sqrt ( (x1-x2)**2 + (y1-y2)**2 ) d13 = sqrt ( (x1-x3)**2 + (y1-y3)**2 ) d32 = sqrt ( (x3-x2)**2 + (y3-y2)**2 ) C x112 = x1 + ( x2 - x1 ) * r1 / d12 x212 = x2 + ( x1 - x2 ) * r2 / d12 y112 = y1 + ( y2 - y1 ) * r1 / d12 y212 = y2 + ( y1 - y2 ) * r2 / d12 x12 = 0.5 * ( x112 + x212 ) y12 = 0.5 * ( y112 + y212 ) C x113 = x1 + ( x3 - x1 ) * r1 / d13 x313 = x3 + ( x1 - x3 ) * r3 / d13 y113 = y1 + ( y3 - y1 ) * r1 / d13 y313 = y3 + ( y1 - y3 ) * r3 / d13 x13 = 0.5 * ( x113 + x313 ) y13 = 0.5 * ( y113 + y313 ) C x332 = x3 + ( x2 - x3 ) * r3 / d32 x232 = x2 + ( x3 - x2 ) * r2 / d32 y332 = y3 + ( y2 - y3 ) * r3 / d32 y232 = y2 + ( y3 - y2 ) * r2 / d32 x32 = 0.5 * ( x332 + x232 ) y32 = 0.5 * ( y332 + y232 ) P(1) = ( x12 + x13 + x32 ) / 3.0 P(2) = ( y12 + y13 + y32 ) / 3.0 ier = 0 C call fminsi ( 2, P, dist3, EPS, FMIN, IER ) C Diagnostic output C write (*,*) ier, FMIN xc = P(1) yc = P(2) C C Radius RC1 = sqrt ( (xc - x1)**2 + (yc - y1)**2 ) - r1 RC2 = sqrt ( (xc - x2)**2 + (yc - y2)**2 ) - r2 RC3 = sqrt ( (xc - x3)**2 + (yc - y3)**2 ) - r3 C write (*,*) rc1, rc2, rc3 RC = ( RC1 + RC2 + RC3 ) / 3.0 end C------------------------------------------- real function dist3 ( n, p ) C Distance from 3 circles integer n real p(n) real x1, x2, x3, y1, y2, y3, r1, r2, r3 common /coor3c/ x1, x2, x3, y1, y2, y3, r1, r2, r3 C x = P(1) y = P(2) C d1 = sqrt (( x-x1)**2 + (y-y1)**2) d2 = sqrt (( x-x2)**2 + (y-y2)**2) d3 = sqrt (( x-x3)**2 + (y-y3)**2) C d1 = d1 - r1 d2 = d2 - r2 d3 = d3 - r3 C dist3 = (d1-d2)**2 + ( d1-d3)**2 + (d2-d3)**2 C & + D1**2 + D2**2 + D3**2 C Diagnostic output C write (*,*) ' DIST3:', x, y, dist3 C end C-------------------------------------- subroutine C20CEN ( M, X, Y, R, I1, I2, & XC1, YC1, RC1, XC2, YC2, RC2 ) C determine the centers and radii of the two circles that touch a C container circle and two other circles that are located inside C the container circle real x(M), Y(M), R(M) real P(2), EPS(2) real x1, x2, x3, y1, y2, y3, r1, r2, r3 common /coor3c/ x1, x2, x3, y1, y2, y3, r1, r2, r3 external dist20 data eps / 0.0, 0.0 / x1 = x(i1) x2 = x(i2) y1 = y(i1) y2 = y(i2) r1 = r(i1) r2 = r(i2) C C Distance of the two interior circles C d12 = sqrt ( (x1-x2)**2 + (y1-y2)**2 ) C Point in the middle of the gap between the two interior circles C DX = X2 - X1 C DY = Y2 - Y1 C x112 = x1 + DX * r1 / d12 C x212 = x2 - DX * r2 / d12 C y112 = y1 + DY * r1 / d12 C y212 = y2 - DY * r2 / d12 C x12 = 0.5 * ( x112 + x212 ) C y12 = 0.5 * ( y112 + y212 ) C construct a normal on the straight line joining the centers of the C two interior circles passing through the gap center and determine C the two intersections with the container circle R0= 1.0 C xc1=( - ((r1 + r2)*(r1 - r2)*(x1 - x2) - (x1**3 - x1**2* & x2 + x2**3) + (sqrt((2*(y1**2 - y1*y2 + y2**2 & + x2**2)*x2 - (x1 - 2*x2)*x1**2 - ( & 2*y1**2 - 2*y1*y2 + y2**2 + 2*x2**2) & *x1)*x1 - (r1**2 - 2*r2**2 - 2*(y1 - y2 & )*y1 - 2*(x1 - x2)*x1)*r1**2 + 2*(( & y1 - y2)**2 + x2**2 + (x1 - 2*x2)*x1 & )*r0**2 - (2*(y1 - y2)*y2 + r2**2 + 2*( & x1 - x2)*x2)*r2**2 - (y1**2 - 2*y1* & y2 + 2*y2**2)*x2**2 - (y1**2 + y2**2)*( & y1 - y2)**2 - x2**4) + x2*y2)*(y1 - y2) & ) + ((y1 - y2)*y1 - x2**2)*x1)/(2*((y1 - y2)**2 & + x2**2) + 2*(x1 - 2*x2)*x1) C yc1=(((y1 - y2)**2 - x1*x2)*(y1 + y2) + x1**2*y1 + x2**2* & y2 + sqrt((2*(y1**2 - y1*y2 + y2**2 + x2**2)*x2 - ( & x1 - 2*x2)*x1**2 - (2*y1**2 - 2*y1*y2 & + y2**2 + 2*x2**2)*x1)*x1 - (r1**2 - 2 & *r2**2 - 2*(y1 - y2)*y1 - 2*(x1 - x2)* & x1)*r1**2 + 2*((y1 - y2)**2 + x2**2 + ( & x1 - 2*x2)*x1)*r0**2 - (2*(y1 - y2)* & y2 + r2**2 + 2*(x1 - x2)*x2)*r2**2 - ( & y1**2 - 2*y1*y2 + 2*y2**2)*x2**2 - (y1 & **2 + y2**2)*(y1 - y2)**2 - x2**4)*(x1 & - x2) - (r1 + r2)*(r1 - r2)*(y1 - y2))/(2*((y1 - & y2)**2 + x2**2) + 2*(x1 - 2*x2)*x1) C C solve minimization problem P(1) = XC1 P(2) = YC1 IER = 0 call fminsi ( 2, P, dist20, eps, fmin, ier ) XC1 = P(1) YC1 = P(2) C Diagnostic output C write (*,*) XC1, YC1, FMIN C RC1 = 1.0 - sqrt ( XC1**2 + Yc1**2 ) C xc2=( - ((r1 + r2)*(r1 - r2)*(x1 - x2) - (x1**3 - x1**2* & x2 + x2**3) - (sqrt((2*(y1**2 - y1*y2 + y2**2 & + x2**2)*x2 - (x1 - 2*x2)*x1**2 - ( & 2*y1**2 - 2*y1*y2 + y2**2 + 2*x2**2) & *x1)*x1 - (r1**2 - 2*r2**2 - 2*(y1 - y2 & )*y1 - 2*(x1 - x2)*x1)*r1**2 + 2*(( & y1 - y2)**2 + x2**2 + (x1 - 2*x2)*x1 & )*r0**2 - (2*(y1 - y2)*y2 + r2**2 + 2*( & x1 - x2)*x2)*r2**2 - (y1**2 - 2*y1* & y2 + 2*y2**2)*x2**2 - (y1**2 + y2**2)*( & y1 - y2)**2 - x2**4) - x2*y2)*(y1 - y2) & ) + ((y1 - y2)*y1 - x2**2)*x1)/(2*((y1 - y2)**2 & + x2**2) + 2*(x1 - 2*x2)*x1) C yc2=(((y1 - y2)**2 - x1*x2)*(y1 + y2) + x1**2*y1 + x2**2* & y2 - sqrt((2*(y1**2 - y1*y2 + y2**2 + x2**2)*x2 - ( & x1 - 2*x2)*x1**2 - (2*y1**2 - 2*y1*y2 & + y2**2 + 2*x2**2)*x1)*x1 - (r1**2 - 2 & *r2**2 - 2*(y1 - y2)*y1 - 2*(x1 - x2)* & x1)*r1**2 + 2*((y1 - y2)**2 + x2**2 + ( & x1 - 2*x2)*x1)*r0**2 - (2*(y1 - y2)* & y2 + r2**2 + 2*(x1 - x2)*x2)*r2**2 - ( & y1**2 - 2*y1*y2 + 2*y2**2)*x2**2 - (y1 & **2 + y2**2)*(y1 - y2)**2 - x2**4)*(x1 & - x2) - (r1 + r2)*(r1 - r2)*(y1 - y2))/(2*((y1 - & y2)**2 + x2**2) + 2*(x1 - 2*x2)*x1) C P(1) = XC2 P(2) = YC2 IER = 0 call fminsi ( 2, P, dist20, eps, fmin, ier ) XC2 = P(1) YC2 = P(2) C Diagnostic output C write (*,*) XC2, YC2, FMIN C RC2 = 1.0 - sqrt ( XC2**2 + Yc2**2 ) C END C--------------------------------------- real function dist20 ( n, p ) C Distance from 2 circles and an outer container circle integer n real p(n) real x1, x2, x3, y1, y2, y3, r1, r2, r3 common /coor3c/ x1, x2, x3, y1, y2, y3, r1, r2, r3 C x = P(1) y = P(2) C d1 = sqrt (( x-x1)**2 + (y-y1)**2) d2 = sqrt (( x-x2)**2 + (y-y2)**2) d0 = 1.0 - sqrt ( x**2 + y**2 ) C d1 = d1 - r1 d2 = d2 - r2 d3 = d0 C dist20 = (d1-d2)**2 + ( d1-d3)**2 + (d2-d3)**2 C Diagnostic output C write (*,*) ' DIST3:', x, y, dist3 C end