# -*- coding: latin-1 -*-

########################################################################################################################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
# Import der benötigten Module
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
########################################################################################################################

# Float bei Division
from __future__ import division

# ABAQUS Module
from abaqusConstants import *
from abaqus import * 
from odbAccess import *
from visualization import *
from caeModules import *
import part
import assembly
import mesh
import regionToolset

import math
import numpy as np

# Parameter für Model
model = 'Test'
part = 'Turm'

# Parameter Kurve
a = 84
b = 209.66

# Funktion für Shell curve
# wir nehmen y=0 an
def shellCurve(z):
	x = math.sqrt( (1 + (z**2)/b**2) * a**2 )
	return x

# Modell erstellen
mdb.Model(name = model)

# Skizze
s1 = mdb.models[model].ConstrainedSketch(name='__profile__', sheetSize=200.0)
geo = s1.geometry
s1.setPrimaryObject(option=STANDALONE)
s1.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 300.0))
s1.FixedConstraint(entity=geo[2])

# Punkte definieren
# hier können beliebige Z-werte definiert werden
zmin = -270
zmax = 60
znum = 100
zlist = np.linspace(zmin,zmax,znum)
xlist = []

# über Z-Punkte iterieren und die zugehörigen X-Werte ausrechnen
for z in zlist:
	xlist.append(shellCurve(z))

# (znum-1) Linien zeichnen
for i in range(znum-1):
	s1.Line( point1=(xlist[i],zlist[i]) , point2=(xlist[i+1],zlist[i+1]) )

# Part erstellen
p = mdb.models[model].Part(name=part, dimensionality=THREE_D, type=DEFORMABLE_BODY)

# Skizze rotieren
p.BaseShellRevolve(sketch=s1, angle=360.0, flipRevolveDirection=OFF)
s1.unsetPrimaryObject()
session.viewports['Viewport: 1'].setValues(displayedObject=p)