Advertisement
AlanSE

Fractal Pattern as SVG for Asteroid Colony

Dec 4th, 2014
466
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ! *********** for pastbin users **********************
  2. ! * Code is supplement to research on space colonies *
  3. ! *   http://gravitationalballoon.blogspot.com/      *
  4. ! ****************************************************
  5. !
  6. ! fractal geometry calcs for rockfill at center of 3 circles
  7. ! going to 6 levels, minimum radius is 0.000445766
  8.  
  9. module level_vals
  10.   double precision :: SA, area, volume
  11.   integer :: iters
  12.   integer :: thecount, gmax
  13.   integer, parameter :: canvas = 500
  14.   double precision, parameter :: pi = 3.1415926535
  15.   double precision :: cutoff = 0.001
  16. end module level_vals
  17.  
  18. program main_prog
  19.   use level_vals
  20.   implicit none  ! main program just calculates power of 3 and passes
  21.                  !  them b/c must be taken as parameter
  22.   integer :: i,j
  23.   double precision, dimension(3,2) :: corners
  24.   character(len=15) :: filename
  25.   character(len=1), dimension(6) :: digit
  26.   double precision, dimension(3) :: radii
  27.  
  28.   ! start by building equalatteral triangle
  29.   corners(1,:) = 0.
  30.   corners(2,:) = (/ cos(60.*pi/180.), sin(60.*pi/180.) /)
  31.   corners(3,:) = (/ 1.0, 0. /)
  32.   radii = 0.5d0
  33.  
  34.   ! open file to dump raw data
  35.   open( unit=11, file="circles.txt", status="replace" )
  36.   write(11,*) ' recursion_level    x_value     y_value     radius '
  37.  
  38.   ! write header markup to the svg output file
  39.   open( unit=31, file='svg/output.svg', status="replace" )
  40.   write(31,'(a)') '<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
  41.   write(31,'(a,I0,a,I0,a)') '<svg xmlns="http://www.w3.org/2000/svg" width="'&
  42.   &,canvas,'" height="',canvas,'">'
  43.  
  44.   ! initialize tallies and call the recursive routine
  45.   SA = 0.
  46.   area = 0.
  47.   volume = 0.
  48.   iters = 0
  49.   thecount = 0
  50.   gmax = 0
  51.   call nook(corners,radii,0)
  52.  
  53.   ! write footer marketup to svg
  54.   write(31,'(a)') "</svg>"
  55.   do i = 0,5
  56.     close(10+i)
  57.   end do
  58.   close(31)
  59.  
  60.   ! summary output of global tallies
  61.   write(*,*) ' '
  62.   write(*,*) '  -- calculated totals --'
  63.   write(*,*) ' '
  64.   write(*,*) '  area    volume   number '
  65.   write(*,*) area, volume, thecount
  66.   write(*,*) ' '
  67.   write(*,*) ' total iterations = ',iters
  68.   write(*,*) ' total rocks      = ',thecount
  69.   write(*,*) ' max recursion    = ',gmax
  70.   write(*,*) ' '
  71.  
  72. end program main_prog
  73.  
  74. recursive subroutine nook(xin,r,g)
  75.   use level_vals
  76.   implicit none
  77.   double precision, dimension(3,2), intent(in) :: xin
  78.   double precision, dimension(3), intent(in) :: r
  79.   integer, intent(in) :: g  ! iteration index
  80.   double precision, dimension(2) :: x_new
  81.   double precision, dimension(3,2) :: xout
  82.   double precision, dimension(3) :: rout
  83.   double precision :: r_new
  84.   integer :: theunit
  85.   interface
  86.     subroutine middle_loc(x_mat,r,x_new,r_new)
  87.       implicit none
  88.       double precision, dimension(3,2), intent(in) :: x_mat
  89.       double precision, dimension(3), intent(in) :: r
  90.       double precision, dimension(2), intent(out) :: x_new
  91.       double precision, intent(out) :: r_new
  92.     end subroutine
  93.   end interface
  94.   double precision :: factor, r_aux(3)
  95.  
  96.   call middle_loc(xin,r,x_new,r_new)
  97.  
  98.   factor = 1.2
  99.   r_new = factor*r_new
  100. !   r_aux = factor*r
  101.  
  102.   if (g==0) then
  103.     write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(1,1)))+1,&
  104.     &'" cy="',floor(wind_y(xin(1,2)))+1,&
  105.     &'" r="',floor(canvas*r(1)/0.5)-1,&
  106.     &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
  107.     write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(2,1))),&
  108.     &'" cy="',floor(wind_y(xin(2,2)))-1,&
  109.     &'" r="',floor(canvas*r(2)/0.5)-1,&
  110.     &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
  111.     write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(xin(3,1)))-1,&
  112.     &'" cy="',floor(wind_y(xin(3,2)))+1,&
  113.     &'" r="',floor(canvas*r(3)/0.5)-1,&
  114.     &'" fill="#aaa" stroke="#000" stroke-width="1" opacity="0.5"/>'
  115.  
  116.   end if
  117.  
  118.  
  119.   write(11,*) g, x_new, r_new
  120.   write(31,'(a,I0,a,I0,a,I0,a)') '<circle cx="',floor(wind_x(x_new(1))),&
  121.   &'" cy="',floor(wind_y(x_new(2))),&
  122.   &'" r="',floor(canvas*r_new/0.5),&
  123.   &'" fill="#29f" stroke="#000" stroke-width="1" opacity="0.5"/>'
  124.  
  125.   r_new = r_new/factor
  126.  
  127.   SA = SA + 4.*pi*(r_new**2)
  128.   area = area + pi*(r_new**2)
  129.   volume = volume + 4.*pi*(r_new**3)/3.
  130.   thecount = thecount + 1
  131.   if (g>gmax) gmax = g
  132.  
  133.   if (r_new>cutoff) then
  134.     xout = xin
  135.     xout(3,:) = x_new(:)
  136.     rout = (/ r(1),r(2),r_new /)
  137.     call nook(xout,rout,g+1)
  138.     xout = xin
  139.     xout(2,:) = x_new(:)
  140.     rout = (/ r(1),r_new,r(3) /)
  141.     call nook(xout,rout,g+1)
  142.     xout = xin
  143.     xout(1,:) = x_new(:)
  144.     rout = (/ r_new,r(2),r(3) /)
  145.     call nook(xout,rout,g+1)
  146.   end if
  147.  
  148.   contains
  149.     double precision function wind_x(x)
  150.       implicit none
  151.       double precision, intent(in) :: x
  152.       wind_x = canvas*(x-0.25)/0.5
  153.     end function wind_x
  154.     double precision function wind_y(y)
  155.       implicit none
  156.       double precision, intent(in) :: y
  157.       wind_y = canvas*y / 0.5
  158.     end function wind_y
  159.  
  160. end subroutine nook
  161.  
  162. subroutine middle_loc(x_mat,r,x_new,r_new)
  163.   use level_vals
  164.   implicit none
  165.   double precision, dimension(3,2), intent(in) :: x_mat
  166.   double precision, dimension(3), intent(in) :: r
  167.   double precision, dimension(2), intent(out) :: x_new
  168.   double precision, intent(out) :: r_new
  169.   double precision, dimension(2) :: xx, move
  170.   double precision, dimension(3) :: radii  ! radius of circle if touching each other circle
  171.   double precision, dimension(3,2) :: x_mid
  172.   double precision :: r_working
  173.   double precision :: dist, v_len
  174.   integer :: i,j,k
  175.  
  176.   ! begin with guess
  177.   do i = 1,3
  178.     j = i+1
  179.     if (j==4) j = 1
  180.     x_mid(j,:) = x_mat(i,:) + &
  181.     &(x_mat(i,:)-x_mat(j,:))/dist(x_mat(i,:),x_mat(j,:))*r(j)
  182.   end do
  183.   xx = (x_mid(1,:)+x_mid(2,:)+x_mid(3,:))/3.0
  184. !   write(*,*) ' start ',xx
  185.  
  186.   k = 0
  187.   do
  188.     ! calculate radii given the 3 proximities
  189.     do j = 1,3
  190.       radii(j) = dist(xx,x_mat(j,:)) - r(j)
  191.     end do
  192.    
  193.     ! identify closest one
  194.     i = minloc(radii,dim=1)
  195.     r_working = radii(i)
  196.    
  197.     ! movement vectors
  198.     move = 0.
  199.     do j = 1,3
  200.       if (j/=i) then
  201.         move = move + (radii(j)-r_working)*(x_mat(j,:)-xx)/(dist(xx,x_mat(j,:))*2.)
  202.       end if
  203.     end do
  204.     xx = xx + move
  205. !     write(*,*) xx,' moved ',v_len(move)
  206. !     write(*,*) '        terms ',v_len(move),r_working
  207.     if (v_len(move)/r_working < 1.0e-5) then
  208.       exit
  209.     end if
  210.     k = k+1
  211.     if (k>40) stop 'failed to converge'
  212.     if (r_working<0) stop 'out of range'
  213.   end do
  214.  
  215.   x_new = xx
  216.   r_new = r_working
  217.   iters = iters + k
  218.  
  219. end subroutine middle_loc
  220.  
  221. double precision function dist(x1,x2)
  222.   implicit none
  223.   double precision, dimension(2), intent(in) :: x1,x2
  224.   dist = sqrt((x1(1)-x2(1))**2+(x1(2)-x2(2))**2)
  225. end function dist
  226.  
  227. double precision function v_len(x1)
  228.   implicit none
  229.   double precision, dimension(2), intent(in) :: x1
  230.   v_len = sqrt((x1(1))**2+(x1(2))**2)
  231. end function v_len
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement