# ----------------------------------------------------------------------------
#
# geometry.py
#
# ----------------------------------------------------------------------------
# latest revision:  - Mar 2020
#
# purpose:          - generates geometry file for full fiber (core + cladding)
#                     with 4 prism and 12 hexa elements
#                     (4 prisms + 4 hexas in core, 8 hexas in cladding).
#
# ----------------------------------------------------------------------------

import math

# Set the length of the fiber (z-direction)
# and open geometry file
LZ =   1.2; f = open("fiber_prism_1_2","w+")
#LZ =   2.4; f = open("fiber_prism_2_4","w+")
#LZ =   4.8; f = open("fiber_prism_4_8","w+")
#LZ =   9.6; f = open("fiber_prism_9_6","w+")
#LZ =  19.2; f = open("fiber_prism_19_2","w+")
#LZ =  38.4; f = open("fiber_prism_38_4","w+")
#LZ =  76.8; f = open("fiber_prism_76_8","w+")
#LZ = 153.6; f = open("fiber_prism_153_6","w+")
#LZ = 307.2; f = open("fiber_prism_307_2","w+")
#LZ = 614.4; f = open("fiber_prism_614_4","w+")
#LZ = 1228.8; f = open("fiber_prism_1228_8","w+")

# Core and cladding size parameters
core = 0.9
clad = 9.0
r_core = core*math.sqrt(2)
r_clad = clad*math.sqrt(2)
print("Fiber length: %.1f" % LZ)
print("Core radius : %.5f" % r_core)
print("Clad radius : %.5f" % r_clad)
inner = core/4
outer = 3.0*core

# Dimension
f.write("3 3 NDIM,MANDIM\n")
f.write("\n")
# Surfaces
f.write("0 NRSURFS\n")
f.write("\n")
# Elements
f.write("4 NRDOMAIN\n")
f.write("\n")
# Points
f.write("34 NRPOINT\n")
f.write("\n")
f.write("Regular             point 1\n")
f.write("-%fd0 -%fd0 0.0d0\n" % (inner, inner))
f.write("\n")
f.write("Regular             point 2\n")
f.write("%fd0 -%fd0 0.0d0\n" % (inner, inner))
f.write("\n")
f.write("Regular             point 3\n")
f.write("%fd0 %fd0 0.0d0\n" % (inner, inner))
f.write("\n")
f.write("Regular             point 4\n")
f.write("-%fd0 %fd0 0.0d0\n" % (inner, inner))
f.write("\n")
f.write("Regular             point 5\n")
f.write("-%fd0 -%fd0 0.0d0\n" % (core, core))
f.write("\n")
f.write("Regular             point 6\n")
f.write("%fd0 -%fd0 0.0d0\n" % (core, core))
f.write("\n")
f.write("Regular             point 7\n")
f.write("%fd0 %fd0 0.0d0\n" % (core, core))
f.write("\n")
f.write("Regular             point 8\n")
f.write("-%fd0 %fd0 0.0d0\n" % (core, core))
f.write("\n")
f.write("Regular             point 9\n")
f.write("-%fd0 -%fd0 0.0d0\n" % (outer, outer))
f.write("\n")
f.write("Regular             point 10\n")
f.write("%fd0 -%fd0 0.0d0\n" % (outer, outer))
f.write("\n")
f.write("Regular             point 11\n")
f.write("%fd0 %fd0 0.0d0\n" % (outer, outer))
f.write("\n")
f.write("Regular             point 12\n")
f.write("-%fd0 %fd0 0.0d0\n" % (outer, outer))
f.write("\n")
f.write("Regular             point 13\n")
f.write("-%fd0 -%fd0 0.0d0\n" % (clad, clad))
f.write("\n")
f.write("Regular             point 14\n")
f.write("%fd0 -%fd0 0.0d0\n" % (clad, clad))
f.write("\n")
f.write("Regular             point 15\n")
f.write("%fd0 %fd0 0.0d0\n" % (clad, clad))
f.write("\n")
f.write("Regular             point 16\n")
f.write("-%fd0 %fd0 0.0d0\n" % (clad, clad))
f.write("\n")
f.write("Regular             point 17\n")
f.write("-%fd0 -%fd0 %.1fd0\n" % (inner, inner, LZ))
f.write("\n")
f.write("Regular             point 18\n")
f.write("%fd0 -%fd0 %.1fd0\n" % (inner, inner, LZ))
f.write("\n")
f.write("Regular             point 19\n")
f.write("%fd0 %fd0 %.1fd0\n" % (inner, inner, LZ))
f.write("\n")
f.write("Regular             point 20\n")
f.write("-%fd0 %fd0 %.1fd0\n" % (inner, inner, LZ))
f.write("\n")
f.write("Regular             point 21\n")
f.write("-%fd0 -%fd0 %.1fd0\n" % (core, core, LZ))
f.write("\n")
f.write("Regular             point 22\n")
f.write("%fd0 -%fd0 %.1fd0\n" % (core, core, LZ))
f.write("\n")
f.write("Regular             point 23\n")
f.write("%fd0 %fd0 %.1fd0\n" % (core, core, LZ))
f.write("\n")
f.write("Regular             point 24\n")
f.write("-%fd0 %fd0 %.1fd0\n" % (core, core, LZ))
f.write("\n")
f.write("Regular             point 25\n")
f.write("-%fd0 -%fd0 %.1fd0\n" % (outer, outer, LZ))
f.write("\n")
f.write("Regular             point 26\n")
f.write("%fd0 -%fd0 %.1fd0\n" % (outer, outer, LZ))
f.write("\n")
f.write("Regular             point 27\n")
f.write("%fd0 %fd0 %.1fd0\n" % (outer, outer, LZ))
f.write("\n")
f.write("Regular             point 28\n")
f.write("-%fd0 %fd0 %.1fd0\n" % (outer, outer, LZ))
f.write("\n")
f.write("Regular             point 29\n")
f.write("-%fd0 -%fd0 %.1fd0\n" % (clad, clad, LZ))
f.write("\n")
f.write("Regular             point 30\n")
f.write("%fd0 -%fd0 %.1fd0\n" % (clad, clad, LZ))
f.write("\n")
f.write("Regular             point 31\n")
f.write("%fd0 %fd0 %.1fd0\n" % (clad, clad, LZ))
f.write("\n")
f.write("Regular             point 32\n")
f.write("-%fd0 %fd0 %.1fd0\n" % (clad, clad, LZ))
f.write("\n")
f.write("Regular             point 33\n")
f.write("0.0d0 0.0d0 0.0d0\n")
f.write("\n")
f.write("Regular             point 34\n")
f.write("0.0d0 0.0d0 %.1fd0\n" % (LZ))
f.write("\n")
# Curves
f.write("81 NRCURVE\n")
f.write("\n")
f.write("QuaCir              curve 1\n")
f.write("1 2\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 2\n")
f.write("2 3\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 3\n")
f.write("4 3\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 4\n")
f.write("1 4\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("Seglin              curve 5\n")
f.write("5 1\n")
f.write("\n")
f.write("Seglin              curve 6\n")
f.write("6 2\n")
f.write("\n")
f.write("Seglin              curve 7\n")
f.write("7 3\n")
f.write("\n")
f.write("Seglin              curve 8\n")
f.write("8 4\n")
f.write("\n")
f.write("QuaCir              curve 9\n")
f.write("5 6\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 10\n")
f.write("6 7\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 11\n")
f.write("8 7\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 12\n")
f.write("5 8\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("Seglin              curve 13\n")
f.write("9 5\n")
f.write("\n")
f.write("Seglin              curve 14\n")
f.write("10 6\n")
f.write("\n")
f.write("Seglin              curve 15\n")
f.write("11 7\n")
f.write("\n")
f.write("Seglin              curve 16\n")
f.write("12 8\n")
f.write("\n")
f.write("QuaCir              curve 17\n")
f.write("9 10\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 18\n")
f.write("10 11\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 19\n")
f.write("12 11\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 20\n")
f.write("9 12\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("Seglin              curve 21\n")
f.write("13 9\n")
f.write("\n")
f.write("Seglin              curve 22\n")
f.write("14 10\n")
f.write("\n")
f.write("Seglin              curve 23\n")
f.write("15 11\n")
f.write("\n")
f.write("Seglin              curve 24\n")
f.write("16 12\n")
f.write("\n")
f.write("QuaCir              curve 25\n")
f.write("13 14\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 26\n")
f.write("14 15\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 27\n")
f.write("16 15\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 28\n")
f.write("13 16\n")
f.write(".0d0, .0d0, .0d0\n")
f.write("\n")
f.write("QuaCir              curve 29\n")
f.write("17 18\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 30\n")
f.write("18 19\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 31\n")
f.write("20 19\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 32\n")
f.write("17 20\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("Seglin              curve 33\n")
f.write("21 17\n")
f.write("\n")
f.write("Seglin              curve 34\n")
f.write("22 18\n")
f.write("\n")
f.write("Seglin              curve 35\n")
f.write("23 19\n")
f.write("\n")
f.write("Seglin              curve 36\n")
f.write("24 20\n")
f.write("\n")
f.write("QuaCir              curve 37\n")
f.write("21 22\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 38\n")
f.write("22 23\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 39\n")
f.write("24 23\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 40\n")
f.write("21 24\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("Seglin              curve 41\n")
f.write("25 21\n")
f.write("\n")
f.write("Seglin              curve 42\n")
f.write("26 22\n")
f.write("\n")
f.write("Seglin              curve 43\n")
f.write("27 23\n")
f.write("\n")
f.write("Seglin              curve 44\n")
f.write("28 24\n")
f.write("\n")
f.write("QuaCir              curve 45\n")
f.write("25 26\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 46\n")
f.write("26 27\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 47\n")
f.write("28 27\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 48\n")
f.write("25 28\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("Seglin              curve 49\n")
f.write("29 25\n")
f.write("\n")
f.write("Seglin              curve 50\n")
f.write("30 26\n")
f.write("\n")
f.write("Seglin              curve 51\n")
f.write("31 27\n")
f.write("\n")
f.write("Seglin              curve 52\n")
f.write("32 28\n")
f.write("\n")
f.write("QuaCir              curve 53\n")
f.write("29 30\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 54\n")
f.write("30 31\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 55\n")
f.write("32 31\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("QuaCir              curve 56\n")
f.write("29 32\n")
f.write(".0d0, .0d0, %.1fd0\n" % (LZ))
f.write("\n")
f.write("Seglin              curve 57\n")
f.write("1 17\n")
f.write("\n")
f.write("Seglin              curve 58\n")
f.write("2 18\n")
f.write("\n")
f.write("Seglin              curve 59\n")
f.write("3 19\n")
f.write("\n")
f.write("Seglin              curve 60\n")
f.write("4 20\n")
f.write("\n")
f.write("Seglin              curve 61\n")
f.write("5 21\n")
f.write("\n")
f.write("Seglin              curve 62\n")
f.write("6 22\n")
f.write("\n")
f.write("Seglin              curve 63\n")
f.write("7 23\n")
f.write("\n")
f.write("Seglin              curve 64\n")
f.write("8 24\n")
f.write("\n")
f.write("Seglin              curve 65\n")
f.write("9 25\n")
f.write("\n")
f.write("Seglin              curve 66\n")
f.write("10 26\n")
f.write("\n")
f.write("Seglin              curve 67\n")
f.write("11 27\n")
f.write("\n")
f.write("Seglin              curve 68\n")
f.write("12 28\n")
f.write("\n")
f.write("Seglin              curve 69\n")
f.write("13 29\n")
f.write("\n")
f.write("Seglin              curve 70\n")
f.write("14 30\n")
f.write("\n")
f.write("Seglin              curve 71\n")
f.write("15 31\n")
f.write("\n")
f.write("Seglin              curve 72\n")
f.write("16 32\n")
f.write("\n")
f.write("Seglin              curve 73\n")
f.write("1 33\n")
f.write("\n")
f.write("Seglin              curve 74\n")
f.write("2 33\n")
f.write("\n")
f.write("Seglin              curve 75\n")
f.write("3 33\n")
f.write("\n")
f.write("Seglin              curve 76\n")
f.write("4 33\n")
f.write("\n")
f.write("Seglin              curve 77\n")
f.write("17 34\n")
f.write("\n")
f.write("Seglin              curve 78\n")
f.write("18 34\n")
f.write("\n")
f.write("Seglin              curve 79\n")
f.write("19 34\n")
f.write("\n")
f.write("Seglin              curve 80\n")
f.write("20 34\n")
f.write("\n")
f.write("Seglin              curve 81\n")
f.write("33 34\n")
f.write("\n")
# Triangles
f.write("8  NRTRIAN\n")
f.write("\n")
f.write("TransTri            triangle 1\n")
f.write("1 2 33\n")
f.write("\n")
f.write("TransTri            triangle 2\n")
f.write("2 3 33\n")
f.write("\n")
f.write("TransTri            triangle 3\n")
f.write("4 3 33\n")
f.write("\n")
f.write("TransTri            triangle 4\n")
f.write("1 4 33\n")
f.write("\n")
f.write("TransTri            triangle 5\n")
f.write("17 18 34\n")
f.write("\n")
f.write("TransTri            triangle 6\n")
f.write("18 19 34\n")
f.write("\n")
f.write("TransTri            triangle 7\n")
f.write("20 19 34\n")
f.write("\n")
f.write("TransTri            triangle 8\n")
f.write("17 20 34\n")
f.write("\n")
# Rectangles
f.write("56 NRRECTA\n")
f.write("\n")
f.write("TraQua              rectangle 1\n")
f.write("5 6 2 1\n")
f.write("\n")
f.write("TraQua              rectangle 2\n")
f.write("6 7 3 2\n")
f.write("\n")
f.write("TraQua              rectangle 3\n")
f.write("8 7 3 4\n")
f.write("\n")
f.write("TraQua              rectangle 4\n")
f.write("5 8 4 1\n")
f.write("\n")
f.write("TraQua              rectangle 5\n")
f.write("9 10 6 5\n")
f.write("\n")
f.write("TraQua              rectangle 6\n")
f.write("10 11 7 6\n")
f.write("\n")
f.write("TraQua              rectangle 7\n")
f.write("12 11 7 8\n")
f.write("\n")
f.write("TraQua              rectangle 8\n")
f.write("9 12 8 5\n")
f.write("\n")
f.write("TraQua              rectangle 9\n")
f.write("13 14 10 9\n")
f.write("\n")
f.write("TraQua              rectangle 10\n")
f.write("14 15 11 10\n")
f.write("\n")
f.write("TraQua              rectangle 11\n")
f.write("16 15 11 12\n")
f.write("\n")
f.write("TraQua              rectangle 12\n")
f.write("13 16 12 9\n")
f.write("\n")
f.write("TraQua              rectangle 13\n")
f.write("21 22 18 17\n")
f.write("\n")
f.write("TraQua              rectangle 14\n")
f.write("22 23 19 18\n")
f.write("\n")
f.write("TraQua              rectangle 15\n")
f.write("24 23 19 20\n")
f.write("\n")
f.write("TraQua              rectangle 16\n")
f.write("21 24 20 17\n")
f.write("\n")
f.write("TraQua              rectangle 17\n")
f.write("25 26 22 21\n")
f.write("\n")
f.write("TraQua              rectangle 18\n")
f.write("26 27 23 22\n")
f.write("\n")
f.write("TraQua              rectangle 19\n")
f.write("28 27 23 24\n")
f.write("\n")
f.write("TraQua              rectangle 20\n")
f.write("25 28 24 21\n")
f.write("\n")
f.write("TraQua              rectangle 21\n")
f.write("29 30 26 25\n")
f.write("\n")
f.write("TraQua              rectangle 22\n")
f.write("30 31 27 26\n")
f.write("\n")
f.write("TraQua              rectangle 23\n")
f.write("32 31 27 28\n")
f.write("\n")
f.write("TraQua              rectangle 24\n")
f.write("29 32 28 25\n")
f.write("\n")
f.write("TraQua              rectangle 25\n")
f.write("1 2 18 17\n")
f.write("\n")
f.write("TraQua              rectangle 26\n")
f.write("2 3 19 18\n")
f.write("\n")
f.write("TraQua              rectangle 27\n")
f.write("4 3 19 20\n")
f.write("\n")
f.write("TraQua              rectangle 28\n")
f.write("1 4 20 17\n")
f.write("\n")
f.write("BilQua              rectangle 29\n")
f.write("5 1 17 21\n")
f.write("\n")
f.write("BilQua              rectangle 30\n")
f.write("6 2 18 22\n")
f.write("\n")
f.write("BilQua              rectangle 31\n")
f.write("7 3 19 23\n")
f.write("\n")
f.write("BilQua              rectangle 32\n")
f.write("8 4 20 24\n")
f.write("\n")
f.write("TraQua              rectangle 33\n")
f.write("5 6 22 21\n")
f.write("\n")
f.write("TraQua              rectangle 34\n")
f.write("6 7 23 22\n")
f.write("\n")
f.write("TraQua              rectangle 35\n")
f.write("8 7 23 24\n")
f.write("\n")
f.write("TraQua              rectangle 36\n")
f.write("5 8 24 21\n")
f.write("\n")
f.write("BilQua              rectangle 37\n")
f.write("9 5 21 25\n")
f.write("\n")
f.write("BilQua              rectangle 38\n")
f.write("10 6 22 26\n")
f.write("\n")
f.write("BilQua              rectangle 39\n")
f.write("11 7 23 27\n")
f.write("\n")
f.write("BilQua              rectangle 40\n")
f.write("12 8 24 28\n")
f.write("\n")
f.write("TraQua              rectangle 41\n")
f.write("9 10 26 25\n")
f.write("\n")
f.write("TraQua              rectangle 42\n")
f.write("10 11 27 26\n")
f.write("\n")
f.write("TraQua              rectangle 43\n")
f.write("12 11 27 28\n")
f.write("\n")
f.write("TraQua              rectangle 44\n")
f.write("9 12 28 25\n")
f.write("\n")
f.write("BilQua              rectangle 45\n")
f.write("13 9 25 29\n")
f.write("\n")
f.write("BilQua              rectangle 46\n")
f.write("14 10 26 30\n")
f.write("\n")
f.write("BilQua              rectangle 47\n")
f.write("15 11 27 31\n")
f.write("\n")
f.write("BilQua              rectangle 48\n")
f.write("16 12 28 32\n")
f.write("\n")
f.write("TraQua              rectangle 49\n")
f.write("13 14 30 29\n")
f.write("\n")
f.write("TraQua              rectangle 50\n")
f.write("14 15 31 30\n")
f.write("\n")
f.write("TraQua              rectangle 51\n")
f.write("16 15 31 32\n")
f.write("\n")
f.write("TraQua              rectangle 52\n")
f.write("13 16 32 29\n")
f.write("\n")
f.write("BilQua              rectangle 53\n")
f.write("1 33 34 17\n")
f.write("\n")
f.write("BilQua              rectangle 54\n")
f.write("2 33 34 18\n")
f.write("\n")
f.write("BilQua              rectangle 55\n")
f.write("3 33 34 19\n")
f.write("\n")
f.write("BilQua              rectangle 56\n")
f.write("4 33 34 20\n")
f.write("\n")
# Prisms
f.write("4 NRPRISM\n")
f.write("\n")
f.write("TIprism             prism 1\n")
f.write("1    1 2 33 17 18 34\n")
f.write("\n")
f.write("TIprism             prism 2\n")
f.write("1    2 3 33 18 19 34\n")
f.write("\n")
f.write("TIprism             prism 3\n")
f.write("1    4 3 33 20 19 34\n")
f.write("\n")
f.write("TIprism             prism 4\n")
f.write("1    1 4 33 17 20 34\n")
f.write("\n")
# Hexahedra
f.write("12 NRHEXAS\n")
f.write("\n")
f.write("TraHex              hexa 1\n")
f.write("2    5 6 2 1 21 22 18 17\n")
f.write("\n")
f.write("TraHex              hexa 2\n")
f.write("2    6 7 3 2 22 23 19 18\n")
f.write("\n")
f.write("TraHex              hexa 3\n")
f.write("2    4 3 7 8 20 19 23 24\n")
f.write("\n")
f.write("TraHex              hexa 4\n")
f.write("2    1 4 8 5 17 20 24 21\n")
f.write("\n")
f.write("TraHex              hexa 5\n")
f.write("3    9 10 6 5 25 26 22 21\n")
f.write("\n")
f.write("TraHex              hexa 6\n")
f.write("3    10 11 7 6 26 27 23 22\n")
f.write("\n")
f.write("TraHex              hexa 7\n")
f.write("3    8 7 11 12 24 23 27 28\n")
f.write("\n")
f.write("TraHex              hexa 8\n")
f.write("3    5 8 12 9 21 24 28 25\n")
f.write("\n")
f.write("TraHex              hexa 9\n")
f.write("4    13 14 10 9 29 30 26 25\n")
f.write("\n")
f.write("TraHex              hexa 10\n")
f.write("4    14 15 11 10 30 31 27 26\n")
f.write("\n")
f.write("TraHex              hexa 11\n")
f.write("4    12 11 15 16 28 27 31 32\n")
f.write("\n")
f.write("TraHex              hexa 12\n")
f.write("4    9 12 16 13 25 28 32 29\n")
f.write("\n")
# Tetrahedra
f.write("0 NRTETRA\n")
f.write("\n")
# Pyramids
f.write("0 NRPYRAM\n")
# Close geometry file
f.close()
