Skip to content

Instantly share code, notes, and snippets.

@siers
Created May 29, 2014
Embed
What would you like to do?
gnuplot-sphere.hs
import Numeric.Matrix hiding (map)
import Data.List.Split
import Text.Printf
type Rotation = Matrix D
rotMatrixX :: D -> Rotation
rotMatrixX rad = fromList . chunksOf 3 $
1 : 0 : 0 :
0 : cos rad : -sin rad :
0 : sin rad : cos rad : []
rotMatrixY :: D -> Rotation
rotMatrixY rad = fromList . chunksOf 3 $
cos rad : 0 : sin rad :
0 : 1 : 0 :
-sin rad : 0 : cos rad : []
rotMatrixZ :: D -> Rotation
rotMatrixZ rad = fromList . chunksOf 3 $
cos rad : -sin rad : 0 :
sin rad : cos rad : 0 :
0 : 0 : 1 : []
rotMatrix :: D -> D -> D -> Rotation
rotMatrix x y z = rotMatrixX x * rotMatrixY y * rotMatrixZ z
type D = Double
type Point = [D]
type Line = [Point]
type Surf = [Line]
circle :: D -> Point
circle x = [cos t, sin t, 0]
where t = x - pi/2
line :: Int -> Point -> D -> Line
line slices (x:y:_) t = [[x, y, -hl], [x, y, hl]]
where l = sin t * 2 * pi / fromIntegral slices
hl = l / 2
slice :: Int -> Surf
slice slices = (line slices =<< circle) `map` [0,0.125..pi]
present :: Line -> String
present ((q:w:e:_):(a:s:d:_):_) = printf "%.3f %.3f %.3f\n%.3f %.3f %.3f\n\n" q w e a s d
rotate :: Rotation -> Surf -> Surf
rotate rotM surf = chunksOf 2 $ rotated
where points = concat surf :: Line
rotated = toList $ fromList points * rotM
ball :: Int -> Surf
ball nSlices = concat $ map (rotate (rotMatrixX (pi/2.0))) >>= (++) $ slices
where rotations = rotMatrixY rads
slices = take nSlices . iterate (rotate rotations) $ slice nSlices
rads = 2 * pi / fromIntegral nSlices
main = mapM_ (putStrLn . present) $ ball 40
#!/usr/bin/env gnuplot
set terminal wxt persist
set mouse
splot "data" using 1:2:3 with lines
@siers
Copy link
Author

siers commented May 29, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment