#!/usr/bin/python

#   Copyright (C) 2005 - Ing. Giorgio Griffon - Venezia

#   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

#disegna una sezione a forma di profilo NACA a 4 cifre suddivisa in quadrilateri e in triangoli importabile nel cgx
#parrebbe funzionare:aumenta i punti intorno al bordo d'attacco
import math
from math import sqrt, atan,pi,sin,cos,pow
#caratteristiche del profilo
corda=0.9#[m]
m=0.02
p=0.4
s=0.09

n=50#n. punti su ogni lato del profilo
np=2*n+1#n. punti del profilo+1

ddee=4#divisioni delle linee esterne
ddii=8#divisioni delle linee interne

nomef='profilo.fbd'

def x(i):#ascissa di un punto della linea media
	g=2#grado della curva:pił alto, maggior addensamento di punti intorno al bordo d'attacco
# 	if i<=n:
# 		x=(i-1.0)/(n-1.0)
# 	else:
# 		x=(2.0*n-i)/(n-1.0)
# 	return x
	if i<=n:
		x=pow(float(i)/n,g)
		print 'sono qui'+str(pow(float(i)/n,g))
	else:
		x=pow(float(2.0*n+1-i)/n,g)
	return x

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

def theta(x):
	if x<=p:
		theta=atan(2*m*(p-x)/(p*p))
	else:
		theta=atan(2*m*(p-x)/(1-p)**2)
	theta=theta+pi/2
	return theta

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

def yc(x,i):#ordinata di un punto del profilo
	if i<=n: yc=(lineamedia(x)+sin(theta(x))*semispessore(x))
	else: yc=(lineamedia(x)-sin(theta(x))*semispessore(x))
	return yc

def xc(x,i):#ascissa di un punto del profilo
	if i<=n: xc=x+cos(theta(x))*semispessore(x)
	else: xc=x-cos(theta(x))*semispessore(x)
	return xc

f = open(nomef,'w')

for i in range(0,np):#generazione dei punti
	xlm=x(i)
	print 'x='+str(xlm)
#	yp1=y(xp1,i)
	xp=corda*xc(xlm,i)
	yp=corda*yc(xlm,i)
	nomep = 'P'+str(i)
	riga = 'PNT '+nomep+'     '+str(xp)+'     '+str(yp)+'     0.00000\n'
	f.write(riga)
	print str(xp)+' - '+str(yp)
nlinee=np-1
for i in range(0,nlinee):#generazione delle linee esterne
	nomel = 'L'+str(i)
	riga = 'LINE '+nomel+' P'+str(i)+' P'+str(i+1)+' 1%(ddee)02d\n' %vars()
	f.write(riga)
riga = 'LINE L'+str(nlinee)+' P'+str(nlinee)+' P0 1%(ddee)02d\n' %vars()#l'ultima linea va dall'ultimo punto al primo
f.write(riga)
for i in range(1,n):#generazione delle linee interne al profilo, cominciano da LI2
	nomel = 'LI'+str(i)
	riga = 'LINE '+nomel+' P'+str(i)+' P'+str(2*n+1-i)+' 1%(ddii)02d\n' %vars()
	f.write(riga)

riga = 'GSUR A0 + BLEND + L0 + LI1 + L'+str(2*n)+'\n'
#prima superficie del naso, triangolare
f.write(riga)
for i in range(1,n-1):#generazione delle superfici, quadrilateri
	riga = 'GSUR A'+str(i)+' + BLEND + L'+str(i)+' + LI'+str(i+1)+' + L'+str(2*n-i)+' - LI'+str(i)+'\n'
	f.write(riga)
riga = 'GSUR A'+str(n-1)+' + BLEND + L'+str(n-1)+' + L'+str(n)+' + L'+str(n+1)+' - LI'+str(n-1)+'\n'
f.write(riga)

f.close()