Thursday, 28 March 2013

Running CADAC4 for the DH-1

This post is the documentation for running the CADAC4 program with the DH1 modules and inputs. All the software has been uploaded to a Google Code project called cadac4-dh1. The initial checkin was the version that I was running from my USB stick under windows,  using the gfortran compiler.

The second checkin was the code that I had under Linux, which can be compiled into an executable using the gfortran compiler.

To compile the main DH1 program, execute 'gfortran CADX3.FOR UTL3.FOR R3_MODULE.FOR'

The output will be an 'a.out' executable which should be renamed as 'DH1'. Running this program then reads the input files and calculates the TRAJ.ASC outputs. Refer to the CADAC4 documentation for full details.

Because I am not certain that compiling the Utility files is working correctly, I have also put the Windows/DOS versions of the utilities (CADIN3.exe etc) into a directory called DOS.

The output files can be plotted in Matlab. Start by opening Matlab in an X session:


First edit the TRAJ.ASC to remove all the headers and column names, t
hen change directory to the DH1 code, and execute the 'PlotGraph.m' script. This will generate a graph similar to the one shown.

The blue line is the altitude which tops out around 250,000 feet. The green lines and the right hand scale show the weight of the vehicle in kg, which starts off at 98,000kg GLOW.

Thursday, 21 March 2013

PyDome Tutorial

Running PyDome

  1.  Install Python. The code was developed under Python 2.7. Packages for all operating systems (Windows/Mac/Linux) are available from 
  2. Create a directory under your users home directory and copy all the PyDome files there. This will be where you run the script from and save the output to.
  3. Update the '' script to set the dome diameter and the geodesic frequency.
  4. Run the '' script from the command line. The default output is to print the results on the screen, so to save it to a text file just redirect the output:

    python > output.txt

An average desktop PC should take around 45 seconds to run for a 6V sphere, with lower frequency calculations taking even less time.
2V Dome results in CAD


The output consists of the list of points with Cartesian coordinates:

*     Points                                             *
Pt77 = (-3782.10808464, 129.73646238, 4656.05271517)
Pt246 = (5547.56463791, 1201.67534466, 1944.35155109)
Pt112 = (4854.10196625, 3526.71151375, 0)
Pt402 = (1854.10196625, 4253.25404176, -3804.22606518)
Pt165 = (-571.430583503, -5647.38560167, -1944.35155109)

As well as the list of edges showing both endpoints and their coordinates:

*     Edges                                              *

Point Pt77 Edge List:
Edge187: Pt74 = [ -3000.0, 974.759088699, 5103.90485011] - Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517]
Edge188: Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517] - Pt76 = [ -3917.72834195, 1272.94710279, 4362.45462008]
Edge200: Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517] - Pt75 = [ -2736.75911987, -209.918005711, 5335.36165135]
Edge202: Pt83 = [ -3464.10161514, -1125.55484451, 4767.92683375] - Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517]
Edge203: Pt79 = [ -4618.03398875, 449.02797658, 3804.22606518] - Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517]
Edge204: Pt77 = [ -3782.10808464, 129.73646238, 4656.05271517] - Pt84 = [ -4428.16927497, -759.490479512, 3976.74377899]

Generated code for CAD Programs

It is easy to modify the code to generate output in a format which can be easily used in a CAD program. I have used CATIA for my examples, as it can run Visual Basic macros to do tasks. I simply created a macro which recorded the creation of two points and a line between them for an edge. Then I added the function ‘Get_CATIA_Desc’ to the Coordinates and Edges classes to print out the data with the right text. To generate all the points and edges in a CATIA part I just copy and paste this output into the VB editor and run the script.

Sample VB code to create a point in CATIA:
Set hybridShapePointCoord369 = hybridShapeFactory1.AddNewPointCoord(1958.86417097, -636.473551394, -5635.40172286)
body1.InsertHybridShape hybridShapePointCoord369
part1.InWorkObject = hybridShapePointCoord369

Program Structure

The advantage of using a language like Python is the Object Oriented (OO) nature of the language, which enables the code to reflect the structure. The objects which are modelled are the geodesic sphere (, icosahedron faces (, edges, (, and coordinates/points ( The main GeoSphere class keeps track of a list of all the points and edges.


Each face of the icosahedron has known coordinates, as determined from the equations in Platonic Solids Tutorial. Given these coordinates we can process each of the 20 faces at a time and use vector maths to calculate the delta vectors from a starting point A to the other two corners B and C. The location of each of the intermediate points can then be determined by scaling the delta vectors according to the frequency, and then iterating to get each row and column. 

Icosahedron with one face divided for a 4V geodesic sphere

For example a frequency 4 dome would have each face divided into 4 rows and 4 columns.

The scaled delta vectors would then be calculated as: 
where f is the geodesic frequency.

With one counter for the row number and one for the column number, we can start at point A and use a loop over these counters to determine each point.

for i in range(0,self.freq_n+1): -- Row number
   for j in range( 1, self.freq_n - i + 1): -- Column number
   # Calculate the new coordinates based on the offset from the initial point (A)
As each face is added to the GeoDome class the coordinates from each face are added to the global list. Once all faces have been added then the list can be checked for duplicate points, as any points on the line where two faces meet will be doubled up.

The final step is to project each point onto the sphere. So far all the calculations have been in Cartesian coordinates (x,y,z), however the Coordinates class also keeps the polar coordinates (theta, phi, r). Each point will then have its polar coordinates radius (‘r’) set to the same value while keeping the angles theta and phi the same. Setting the radius will also recalculate the Cartesian coordinates for each point.
Icosahedron with sphere for projecting points

Advanced Modifications


The user parameters are defined in a small section at the beginning of the '' script. Currently there are three parameters, the frequency of the geodesic (frequency_n), the radius of the dome (R_mm), and a flag for sphere or dome calculations. The coordinates should be entered in mm as this gives better readability in the final results.

To calculate the coordinates of a 6 metre (6000mm) radius dome of frequency 6, set these variables to the corresponding values as shown:

R_mm = 6000
frequency_n = 6
Dome_calc = True

There are potentially lots of different methods of dividing the faces of the icosahedron. Traditional geodesic domes or spheres divide each flat face into equal lengths. When the points are projected onto
the sphere they then have different lengths. An alternative method is to divide the arc between the edge vertices into equal angles. In order to change the calculation method for each face,
modify the file in the function Add_Face. This section of the code takes the 3 corner vertices and then calculates all the points within that face. Just replace the function 'Get_Edges_Equal_Distance'
with an equivalent function.

        # Use this function to divide faces by equal distance
        edge_list = F1.Get_Edges_Equal_Distance()

To Do List

For each hub point calculate the angle of at which each edge meets at the hub. This would be both the angle between the edges and the angles of the edge from a flat plane.

Tuesday, 12 March 2013

Geodesic Domes using Python

Over the last few months I have been learning some of the maths behind geodesic domes. I was specifically interested in being able to generate the coordinates for all the points in a way which had enough accuracy to be able to use in a CAD program.

While I was investigating how the vertex points are generated I found this great site with all the information necessary to calculate the main points of the Icosahedron. From there the algorithm is fairly simple to explain, if a little harder to implement! Each of the 20 triangular faces is divided into smaller triangles according to the frequency of the dome. For example a 5V dome would have each face divided into 5 rows of equally sized triangles. From there each point is projected outwards onto a sphere to create the dome arrangement.

To make this easier and faster, I have written a small program in Python to do exactly this. At the moment it takes in the radius of the geodesic sphere and the frequency, and outputs the coordinates and edges. To make it easier I also made the printed output of the coordinates and edges as a piece of text which can be copied into the CATIA VB script tool, and then run as a VB script to automatically create the points and lines.

If you just want to print the cartesian coordinates then the printing of this CATIA code can be commented out.

As an example, here is a picture in CATIA of a frequency 9 dome with a 10 metre diameter.

Here is another one of a smaller frequency 5 sphere with a 4 metre diameter.

In order to share the code I have created a project on Google Code in a Subversion repository. You can download the code from PyDom or checkout directly with Subversion using the command:

svn checkout pydome-read-only

Please let me know any suggestions for improvements!

The biggest enhancement to this script in the future will be the ability to calculate the angles at which the edges meet at each point. This would be crucial to enable the construction of the hubs to anchor all the edges.

Please let me know of any suggestions or fixes, keeping in mind that this is a work in progress!