import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
geompy = geomBuilder.New(salome.myStudy)
gg = salome.ImportComponentGUI("GEOM")
import sys
from math import *
from StringIO import *
# This program makes a turbine profile
r=145 # radius of a rotor
stat_r=155 # radius of a stator
height=1.0 # height of a turbine - 2D case, just one cell 
w=h=2*r
p0 = geompy.MakeVertex(0., 0., 0.)
px = geompy.MakeVertex(r * stat_r, 0., 0.)
py = geompy.MakeVertex(0.,r * stat_r, 0.)
pz = geompy.MakeVertex(0.,0.,height) # height of a turbine
v0z  = geompy.MakeVector(p0, pz)
id_v0z=geompy.addToStudy(v0z,"Vector_z")
gg.createAndDisplayGO(id_v0z)
n=8 # number of vertices in the profile
factor=0.55 #factor to decrease entrance of vortex 
factor_02=0.9 #factor to increase size of internal cavities
angle_01 = factor*2*pi/n
angle_02 = (2*pi/n)*(1-factor)
cosal=cos(angle_01)
sinal=sin(angle_01)
cosal_02=cos(angle_02)
sinal_02=sin(angle_02)
x2 = b = 0
y1=r*cosal
y2=r
x3=r*sinal
y3=r*cosal
points_02 = []
points_03 = []
arcs = []
while b < n*2 :
    a=StringIO()
    points_02.append(geompy.MakeVertex(x2,y2,0))
    a.write('p_02_%i' % b)
    geompy.addToStudy(points_02[b], a.getvalue())
    a.close()    
    a=StringIO()
    points_03.append(geompy.MakeVertex(x3,y3,0))
    a.write('p_03_%i' % b)
    geompy.addToStudy(points_03[b], a.getvalue())
    a.close()
    if (b % 2) == 0:
        arcs.append(geompy.MakeArcCenter(p0, points_02[b], points_03[b], 0))
        x2 = x3 
        y2 = y3 
        tx3=x3
        x3=x3*cosal_02+y3*sinal_02
        y3=-tx3*sinal_02+y3*cosal_02
    else:
        p_center=geompy.MakeVertex(factor_02*x2*cos(angle_02/2)+factor_02*y2*sin(angle_02/2),-factor_02*x2*sin(angle_02/2)+factor_02*y2*cos(angle_02/2),0)
        arcs.append(geompy.MakeArcCenter(p_center, points_02[b], points_03[b], 1))
        x2=x3 
        y2=y3 
        tx3=x3
        x3=x3*cosal+y3*sinal
        y3=-tx3*sinal+y3*cosal
    a=StringIO()
    a.write('arc_%i' % b)
    geompy.addToStudy(arcs[b], a.getvalue())
    a.close() 
    b = b+1
compound_01=geompy.MakeCompound(arcs)
id_compound_01=geompy.addToStudy(compound_01,"Compound_01")
gg.createAndDisplayGO(id_compound_01)
inre_circle=geompy.MakeCircle(p0,v0z,14.5) # radius of internal circle to put turbine on a shaft
id_circle=geompy.addToStudy(inre_circle,"inre_circle")
gg.createAndDisplayGO(id_circle)
face_01=geompy.MakeFaces([inre_circle,compound_01],1)
id_face_01=geompy.addToStudy(face_01,"Face_01")
gg.createAndDisplayGO(id_face_01)
edge = geompy.MakeEdge(p0, pz)
compound_02=geompy.MakeTranslationVector(compound_01,v0z)
id_compound_02=geompy.addToStudy(compound_02,"Compound_02") # internal boundary for upper face
gg.createAndDisplayGO(id_compound_02)
ute_circle_01=geompy.MakeCircle(p0,v0z,165) # diameter of vortex tube 
id_ute_circle_01=geompy.addToStudy(ute_circle_01,"Ute_circle_01")
gg.createAndDisplayGO(id_ute_circle_01)
ute_circle_02=geompy.MakeTranslationVector(ute_circle_01,v0z)
id_ute_circle_02=geompy.addToStudy(ute_circle_02,"Ute_circle_02") # external boundary for upper face
gg.createAndDisplayGO(id_ute_circle_02)
pipe_01=geompy.MakePipe(compound_01,edge)
id_pipe_01=geompy.addToStudy(pipe_01,"Pipe_01") # rotor wall
gg.createAndDisplayGO(id_pipe_01)
pipe_02=geompy.MakePipe(ute_circle_01,edge)
id_pipe_02=geompy.addToStudy(pipe_02,"Pipe_02") # stator wall
gg.createAndDisplayGO(id_pipe_02)
face_02=geompy.MakeFaces([compound_01,ute_circle_01],1)
id_face_02=geompy.addToStudy(face_02,"Face_02") # lower base for volume of air between stator and rotor
gg.createAndDisplayGO(id_face_02)
face_03=geompy.MakeFaces([compound_02,ute_circle_02],1) # upper base for volume of air between stator and rotor
id_face_03=geompy.addToStudy(face_03,"Face_03")
gg.createAndDisplayGO(id_face_03)
v_shell=geompy.MakeShell([face_02,pipe_01,pipe_02,face_03])
id_v_shell=geompy.addToStudy(v_shell,"V_Shell")
gg.createAndDisplayGO(id_v_shell)
v_air=geompy.MakeSolid([v_shell])
id_v_air=geompy.addToStudy(v_air,"Air")
gg.createAndDisplayGO(id_v_air)