AlanSE

Fractal Pattern as SVG for Asteroid Colony

Dec 4th, 2014
228
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
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×