import gmsh
import sys
import numpy as np

# set up parameters (eventually will read from file).
th_l = 10.0
th_u = 40.0

phioff = 160.0
phi = 20.0

infrac = 0.0
outfrac = 0.7
thickness = 1.0

radius = 1.0

# create rectangular coords of point from r, theta, phi.
def polars(rad,th,phi):
    x = rad*np.cos(np.radians(th))*np.cos(np.radians(phi))
    y = rad*np.cos(np.radians(th))*np.sin(np.radians(phi))
    z = rad*np.sin(np.radians(th))
    return x,y,z

# start of main code
gmsh.initialize(sys.argv)
gmsh.model.add("sphere_cut")
sph1 = gmsh.model.occ.addSphere(0,0,0, radius, -1, 0, np.pi/2, 2*np.pi)
sph2 = gmsh.model.occ.addSphere(0,0,0, radius+thickness, -1, 0, np.pi/2, 2*np.pi)

print(type(sph1))
print(sph1,sph2)

shell = gmsh.model.occ.cut([(3,sph2)],[(3,sph1)],5)
gmsh.model.occ.synchronize()

print(type(shell))
print('shell = ',shell)

shellb = gmsh.model.getBoundary([(3,5)])

print(type(shellb))
print('shellb = ',shellb)

epi = shellb[0][1]
top = shellb[1][1]
endo = shellb[2][1]

print("epi: %d, top: %d, endo: %d"%(epi,top,endo))

origin = gmsh.model.occ.addPoint(0.,0.,0.,0,-1)

x,y,z = polars(radius+infrac*thickness,th_l,phioff-phi)
pll = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+infrac*thickness,th_u,phioff-phi)
pul = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+infrac*thickness,th_u,phioff+phi)
pur = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+infrac*thickness,th_l,phioff+phi)
plr = gmsh.model.occ.addPoint(x,y,z,0,-1)

c11 = gmsh.model.occ.addCircleArc(pll,origin,pul,-1)
c21 = gmsh.model.occ.addCircleArc(pul,origin,pur,-1)
c31 = gmsh.model.occ.addCircleArc(pur,origin,plr,-1)
c41 = gmsh.model.occ.addCircleArc(plr,origin,pll,-1)

inwire = gmsh.model.occ.addWire([c11,c21,c31,c41],-1)

x,y,z = polars(radius+outfrac*thickness,th_l,phioff-phi)
pll = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+outfrac*thickness,th_u,phioff-phi)
pul = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+outfrac*thickness,th_u,phioff+phi)
pur = gmsh.model.occ.addPoint(x,y,z,0,-1)
x,y,z = polars(radius+outfrac*thickness,th_l,phioff+phi)
plr = gmsh.model.occ.addPoint(x,y,z,0,-1)

c12 = gmsh.model.occ.addCircleArc(pll,origin,pul,-1)
c22 = gmsh.model.occ.addCircleArc(pul,origin,pur,-1)
c32 = gmsh.model.occ.addCircleArc(pur,origin,plr,-1)
c42 = gmsh.model.occ.addCircleArc(plr,origin,pll,-1)

outwire = gmsh.model.occ.addWire([c12,c22,c32,c42],-1)

in_surf = gmsh.model.occ.addSurfaceFilling(inwire,-1)
out_surf = gmsh.model.occ.addSurfaceFilling(outwire,-1)

isch = gmsh.model.occ.addThruSections([inwire,outwire],-1,1)
gmsh.model.occ.synchronize()
print(type(isch))
print('isch =',isch)
ischb = gmsh.model.getBoundary([(3,isch[0][1])])

print(type(ischb))
print('ischb = ',ischb)

gmsh.model.occ.remove([(2,in_surf),(2,out_surf)])
norm_tiss = gmsh.model.occ.cut(shell[0],isch,-1,0,0)
print('normal: ',norm_tiss)
gmsh.model.occ.synchronize()
isch_tiss = gmsh.model.occ.intersect(shell[0],isch,-1,1,0)
print('ischaemic: ',isch_tiss)

#gmsh.model.removeEntities([(3,sph)])
#gmsh.model.removeEntities([(2,2), (2,4), (2,6)], True)
#gmsh.model.occ.removeAllDuplicates
gmsh.model.occ.synchronize()
gmsh.fltk.run()
gmsh.finalize()