#sintetizza una pala sulla base dei dati di progetto: usa il vortice libero, profili NACA a 4 cifre
#   Copyright (C) 2003 - Giorgio Griffon

#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.

#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.

#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


import Blender
import math
from Blender import NMesh,Object
from math import sqrt,sin,cos,atan,floor,ceil,pi

#####################################################
#DATI DI PROGETTO
p=0.4	#ascissa punto massimo
si=0.09 #semispessore al mozzo
se=0.04	#sem. periferico
ri=0.072	#raggio interno pala (rispetto asse)
re=0.147	#raggio esterno pala (rispetto asse)
Q=0.28 #portata [m^3/s]
ht=5 #prevalenza teorica [m]
ngiri=1450 #giri/min
Z=4 #n. pale
Csi=0.85 #Cs al mozzo, occhio allo stallo
Cse=0.3 #Cs esterno

#DATI DI PROGETTO ACCESSORI
csi=0.85 #coeff. d'ostruzione palare
ns=0.5 #per schiere deceleranti, pompe; per turbine impostare ns=1
xC_rel=0.45 #ascissa centro di svergolamento (1/corda)

#PARAMETRI DELL'APPROSSIMAZIONE GEOMETRICA
dt=0.01	#distanza tra punti successivi dello stesso profilo
l=40 #n. profili accatastati (o va da 0 a l-1)
######################################################

#costanti calcolate, non modificare
g=9.81
deltaz=(re-ri)/(l-1) #[m]distanza tra i piani di giacitura di profili successivi
ght=g*ht
omega=ngiri*2*pi/60
cm=Q/(pi*csi*(re**2-ri**2))

def Cs(r): #variazione di Cs con il raggio, qui lineare
	c=Csi+(Cse-Csi)*(r-ri)/(re-ri)
	return c

def s(r):	#variazione del semispessore col raggio, qui lineare
	sem=si+(se-si)*(r-ri)/(re-ri)
	return sem

def x(t):
	if t<=1: x=corda*t
	else: x=corda*(2-t)
	return x

def lineamedia(t):
	if t<=p: ly=(m/p**2)*(2*p*t-t**2)
	else: ly=(m/(1-p)**2)*(1-2*p+2*p*t-t**2)
	return ly

def semispessore(t):
	deltay=(s(r)/0.2)*(0.2969*sqrt(t)-0.1260*t-0.3516*t**2+0.2843*t**3-0.1015*t**4)
	return deltay

def y(t):
	if t<=1: y=corda*(lineamedia(t)+semispessore(t))
	else: y=corda*(lineamedia(2-t)-semispessore(2-t))
	return y

listax=[]
listay=[]
listaz=[]
listaCu=[]

#scrive le liste di coordinate
print('informazioni');
for o in range(0,l):
	r=ri+deltaz*o
	u=omega*r
	deltaCu=ght/u
	listaCu.append(deltaCu)
	winf=sqrt((u-deltaCu/2)**2+cm**2)
	beta1=atan(cm/u)
	beta2=atan(cm/(u-deltaCu))
	theta=beta2-beta1
	passo=2*pi*r/Z
	corda=2*deltaCu*passo/(winf*Cs(r))
	ms=0.92*(p**2)+(90-beta2*180/pi)/500
	thetac=theta/(1-ms*((passo/corda)**ns))
	m=p*(1-p)*thetac/2 #ordinata massima del profilo
	beta2c=beta1+thetac
	theta1=atan(2*m/p)
	gamma=beta1+theta1

	t=0
	i=0
	xC=xC_rel*corda
	yC=y(xC_rel)
	while t<2:
		xP=x(t)
		yP=y(t)	#coord. profilo piano non svergolato
		xP1=cos(gamma)*(xP-xC)+sin(gamma)*(yP-yC)
		yP1=-sin(gamma)*(xP-xC)+cos(gamma)*(yP-yC)+yC #coord. profilo piano svergolato
		xP2=r*sin(xP1/r)
		zP2=r-r*(1-cos(xP1/r)) #coord. profilo incurvato svergolato
		listax.append(xP2)
		listay.append(yP1)
		listaz.append(zP2)
		i=i+1
		t=i*dt
	o=o+1

n=i #n. punti di ogni profilo
#crea i vertici in base alle coordinate
girante=NMesh.GetRaw()
for j in range (0,l*n):
	nodo=NMesh.Vert(listax[j],listay[j],listaz[j])
	girante.verts.append(nodo)

#crea le facce del solido
for j in range(0,n*(l-1)): #il solido ha n(l-1) facce laterali
	faccia=NMesh.Face()
	if n*int((j+1)/n)<(j+1):
		faccia.v.append(girante.verts[j])
		faccia.v.append(girante.verts[j+n])
		faccia.v.append(girante.verts[j+n+1])
		faccia.v.append(girante.verts[j+1])
	else:
		faccia.v.append(girante.verts[j])
		faccia.v.append(girante.verts[j+n])
		faccia.v.append(girante.verts[j+1])
		faccia.v.append(girante.verts[j-n+1])
	girante.faces.append(faccia)

NMesh.PutRaw(girante,"funzione",1)

#Ricomincia il calcolo per la pala statorica
#####################################################
#####################################################
#DATI DI PROGETTO
pS=0.4	#ascissa punto massimo
siS=0.09 #semispessore al mozzo
seS=0.04	#sem. periferico
riS=0.066
reS=0.147
ZS=5 #n. pale
CsiS=1.0 #Cs al mozzo, occhio stallo
CseS=0.5 #Cs esterno

#DATI DI PROGETTO ACCESSORI
csiS=0.85 #coeff. d'ostruzione palare
xC_relS=0.45 #ascissa centro di svergolamento (1/corda)
distanza=0.3	#distanza assiale a cui collocare le pale statoriche, da aggiustare
######################################################

#costanti calcolate, non modificare
deltaz=(reS-riS)/(l-1) #[m]distanza tra i piani di giacitura di profili successivi

def CsS(r): #variazione di Cs con il raggio, qui lineare
	c=CsiS+(CseS-CsiS)*(r-riS)/(reS-riS)
	return c

def sS(r):	#variazione del semispessore col raggio, qui lineare
	sem=siS+(seS-siS)*(r-riS)/(reS-riS)
	return sem

def xS(t):
	if t<=1: x=corda*t
	else: x=corda*(2-t)
	return x

def lineamediaS(t):
	if t<=pS: ly=(m/pS**2)*(2*pS*t-t**2)
	else: ly=(m/(1-pS)**2)*(1-2*pS+2*pS*t-t**2)
	return ly

def semispessoreS(t):
	deltay=(sS(r)/0.2)*(0.2969*sqrt(t)-0.1260*t-0.3516*t**2+0.2843*t**3-0.1015*t**4)
	return deltay

def yS(t):
	if t<=1: y=corda*(lineamediaS(t)+semispessoreS(t))
	else: y=corda*(lineamediaS(2-t)-semispessoreS(2-t))
	return y

listax=[]
listay=[]
listaz=[]

#scrive le liste di coordinate
for o in range(0,l):
	r=riS+deltaz*o
	Cu=listaCu[o]
	cinf=sqrt((Cu/2)**2+cm**2)
	alfa2=atan(cm/Cu)
	alfa3=pi/2
	theta=alfa3-alfa2
	passo=2*pi*r/ZS
	corda=2*Cu*passo/(cinf*CsS(r))
	ms=0.92*(pS**2)+(90-alfa3*180/pi)/500
	thetac=theta/(1-ms*((passo/corda)**ns))
	m=pS*(1-pS)*thetac/2 #ordinata massima del profilo
	alfa3c=alfa2+thetac
	theta1=atan(2*m/pS)
	gamma=alfa2+theta1

	t=0
	i=0
	xC=xC_relS*corda
	yC=yS(xC_relS)
	while t<2:
		xP=xS(t)
		yP=yS(t)	#coord. profilo piano non svergolato
		xP1=cos(gamma)*(xP-xC)+sin(gamma)*(yP-yC)
		yP1=-sin(gamma)*(xP-xC)+cos(gamma)*(yP-yC)+yC #coord. profilo piano svergolato
		xP2=r*sin(xP1/r)
		zP2=r-r*(1-cos(xP1/r)) #coord. profilo incurvato svergolato
		listax.append(-xP2)
		listay.append(yP1-distanza)
		listaz.append(zP2)
		i=i+1
		t=i*dt
	o=o+1

n=i #n. punti di ogni profilo
#crea i vertici in base alle coordinate
stat=NMesh.GetRaw()
for j in range (0,l*n):
	nodo=NMesh.Vert(listax[j],listay[j],listaz[j])
	stat.verts.append(nodo)

#crea le facce del solido
for j in range(0,n*(l-1)): #il solido ha n(l-1) facce laterali
	faccia=NMesh.Face()
	if n*int((j+1)/n)<(j+1):
		faccia.v.append(stat.verts[j])
		faccia.v.append(stat.verts[j+n])
		faccia.v.append(stat.verts[j+n+1])
		faccia.v.append(stat.verts[j+1])
	else:
		faccia.v.append(stat.verts[j])
		faccia.v.append(stat.verts[j+n])
		faccia.v.append(stat.verts[j+1])
		faccia.v.append(stat.verts[j-n+1])
	stat.faces.append(faccia)

NMesh.PutRaw(stat,"funzione2",1)
Blender.Redraw()