Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- ! *********** for pastbin users **********************
- ! * Code is supplement to research on space colonies *
- ! * http://gravitationalballoon.blogspot.com/ *
- ! ****************************************************
- !
- ! fractal geometry calcs for rockfill at center of 3 circles
- ! going to 6 levels, minimum radius is 0.000445766
- module level_vals
- double precision :: SA, area, volume
- integer :: iters
- integer :: thecount, gmax
- integer, parameter :: canvas = 500
- double precision, parameter :: pi = 3.1415926535
- double precision :: cutoff = 0.001
- end module level_vals
- program main_prog
- use level_vals
- implicit none ! main program just calculates power of 3 and passes
- ! them b/c must be taken as parameter
- integer :: i,j
- double precision, dimension(3,2) :: corners
- character(len=15) :: filename
- character(len=1), dimension(6) :: digit
- double precision, dimension(3) :: radii
- ! start by building equalatteral triangle
- corners(1,:) = 0.
- corners(2,:) = (/ cos(60.*pi/180.), sin(60.*pi/180.) /)
- corners(3,:) = (/ 1.0, 0. /)
- radii = 0.5d0
- ! open file to dump raw data
- open( unit=11, file="circles.txt", status="replace" )
- write(11,*) ' recursion_level x_value y_value radius '
- ! write header markup to the svg output file
- open( unit=31, file='svg/output.svg', status="replace" )
- write(31,'(a)') '<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
- write(31,'(a,I0,a,I0,a)') '<svg xmlns="http://www.w3.org/2000/svg" width="'&
- &,canvas,'" height="',canvas,'">'
- ! initialize tallies and call the recursive routine
- SA = 0.
- area = 0.
- volume = 0.
- iters = 0
- thecount = 0
- gmax = 0
- call nook(corners,radii,0)
- ! write footer marketup to svg
- write(31,'(a)') "</svg>"
- do i = 0,5
- close(10+i)
- end do
- close(31)
- ! summary output of global tallies
- write(*,*) ' '
- write(*,*) ' -- calculated totals --'
- write(*,*) ' '
- write(*,*) ' area volume number '
- write(*,*) area, volume, thecount
- write(*,*) ' '
- write(*,*) ' total iterations = ',iters
- write(*,*) ' total rocks = ',thecount
- write(*,*) ' max recursion = ',gmax
- write(*,*) ' '
- end program main_prog
- recursive subroutine nook(xin,r,g)
- use level_vals
- implicit none
- double precision, dimension(3,2), intent(in) :: xin
- double precision, dimension(3), intent(in) :: r
- integer, intent(in) :: g ! iteration index
- double precision, dimension(2) :: x_new
- double precision, dimension(3,2) :: xout
- double precision, dimension(3) :: rout
- double precision :: r_new
- integer :: theunit
- interface
- subroutine middle_loc(x_mat,r,x_new,r_new)
- implicit none
- double precision, dimension(3,2), intent(in) :: x_mat
- double precision, dimension(3), intent(in) :: r
- double precision, dimension(2), intent(out) :: x_new
- double precision, intent(out) :: r_new
- end subroutine
- end interface
- double precision :: factor, r_aux(3)
- call middle_loc(xin,r,x_new,r_new)
- factor = 1.2
- r_new = factor*r_new
- ! r_aux = factor*r
- if (g==0) then
- write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(1,1)))+1,&
- &'" cy="',floor(wind_y(xin(1,2)))+1,&
- &'" r="',floor(canvas*r(1)/0.5)-1,&
- &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
- write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(2,1))),&
- &'" cy="',floor(wind_y(xin(2,2)))-1,&
- &'" r="',floor(canvas*r(2)/0.5)-1,&
- &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
- write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(3,1)))-1,&
- &'" cy="',floor(wind_y(xin(3,2)))+1,&
- &'" r="',floor(canvas*r(3)/0.5)-1,&
- &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
- end if
- write(11,*) g, x_new, r_new
- write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(x_new(1))),&
- &'" cy="',floor(wind_y(x_new(2))),&
- &'" r="',floor(canvas*r_new/0.5),&
- &'" fill="#29f" stroke="#000" stroke-width="1" opacity="0.5"/>'
- r_new = r_new/factor
- SA = SA + 4.*pi*(r_new**2)
- area = area + pi*(r_new**2)
- volume = volume + 4.*pi*(r_new**3)/3.
- thecount = thecount + 1
- if (g>gmax) gmax = g
- if (r_new>cutoff) then
- xout = xin
- xout(3,:) = x_new(:)
- rout = (/ r(1),r(2),r_new /)
- call nook(xout,rout,g+1)
- xout = xin
- xout(2,:) = x_new(:)
- rout = (/ r(1),r_new,r(3) /)
- call nook(xout,rout,g+1)
- xout = xin
- xout(1,:) = x_new(:)
- rout = (/ r_new,r(2),r(3) /)
- call nook(xout,rout,g+1)
- end if
- contains
- double precision function wind_x(x)
- implicit none
- double precision, intent(in) :: x
- wind_x = canvas*(x-0.25)/0.5
- end function wind_x
- double precision function wind_y(y)
- implicit none
- double precision, intent(in) :: y
- wind_y = canvas*y / 0.5
- end function wind_y
- end subroutine nook
- subroutine middle_loc(x_mat,r,x_new,r_new)
- use level_vals
- implicit none
- double precision, dimension(3,2), intent(in) :: x_mat
- double precision, dimension(3), intent(in) :: r
- double precision, dimension(2), intent(out) :: x_new
- double precision, intent(out) :: r_new
- double precision, dimension(2) :: xx, move
- double precision, dimension(3) :: radii ! radius of circle if touching each other circle
- double precision, dimension(3,2) :: x_mid
- double precision :: r_working
- double precision :: dist, v_len
- integer :: i,j,k
- ! begin with guess
- do i = 1,3
- j = i+1
- if (j==4) j = 1
- x_mid(j,:) = x_mat(i,:) + &
- &(x_mat(i,:)-x_mat(j,:))/dist(x_mat(i,:),x_mat(j,:))*r(j)
- end do
- xx = (x_mid(1,:)+x_mid(2,:)+x_mid(3,:))/3.0
- ! write(*,*) ' start ',xx
- k = 0
- do
- ! calculate radii given the 3 proximities
- do j = 1,3
- radii(j) = dist(xx,x_mat(j,:)) - r(j)
- end do
- ! identify closest one
- i = minloc(radii,dim=1)
- r_working = radii(i)
- ! movement vectors
- move = 0.
- do j = 1,3
- if (j/=i) then
- move = move + (radii(j)-r_working)*(x_mat(j,:)-xx)/(dist(xx,x_mat(j,:))*2.)
- end if
- end do
- xx = xx + move
- ! write(*,*) xx,' moved ',v_len(move)
- ! write(*,*) ' terms ',v_len(move),r_working
- if (v_len(move)/r_working < 1.0e-5) then
- exit
- end if
- k = k+1
- if (k>40) stop 'failed to converge'
- if (r_working<0) stop 'out of range'
- end do
- x_new = xx
- r_new = r_working
- iters = iters + k
- end subroutine middle_loc
- double precision function dist(x1,x2)
- implicit none
- double precision, dimension(2), intent(in) :: x1,x2
- dist = sqrt((x1(1)-x2(1))**2+(x1(2)-x2(2))**2)
- end function dist
- double precision function v_len(x1)
- implicit none
- double precision, dimension(2), intent(in) :: x1
- v_len = sqrt((x1(1))**2+(x1(2))**2)
- end function v_len
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement