"""
python script that generates mesh for Wei et al geometry
"""

import numpy
import nanopores.py4gmsh.basic
import nanopores.py4gmsh.extra
from nanopores.py4gmsh import *
from params_geo import *
from warnings import warn

def get_geo(x0 = None, **params):
    """
    writes a 3d geo file for an axissymmetric geometry for Wei et al
    'Stochastic sensing ...'
    _________
    |        |
    |        |
    |   _____|
    |   \____|
    |        |
    |        |
    |________|   *rotated around z-axis

    """

    reload(nanopores.py4gmsh.basic)
    reload(nanopores.py4gmsh.extra)
    globals().update(params)

    # define additional geo variables
    # mesh generation should also work for no SAM layer
    if lsam < tolc or lsam is None:
        sam = None
    else:
        sam = True

    l0 =  lau +lsin +lsam

    angle2 = angle/2.0
    tan = numpy.tan(angle2*numpy.pi/180)
    sin = numpy.sin(angle2*numpy.pi/180)
    cos = numpy.cos(angle2*numpy.pi/180)

    r1 = r0 + l0*tan
    rsam = r0 + lsam/cos
    rsin = r0 + lsam/cos + rlau
    Rx = R

    X_Fluid_up = numpy.array([
        [0, 0, Rz],
        [Rx, 0, Rz],
        [Rx, 0, l0/2],
    ])

    X_Fluid_low = numpy.array([
        [Rx, 0, -l0/2],
        [Rx, 0, -Rz],
        [0, 0, -Rz],
    ])

    X_Fluid_ctr = numpy.array([
        [0, 0, -l0/2],
        [0, 0, -l0/6],
        [0, 0, +l0/6],
        [0, 0, +l0/2],
    ])

    X_SAM_ctr = numpy.array([
        [r1, 0, l0/2],
        [(2*r1+r0)/3, 0, l0/6],
        [(r1+2*r0)/3, 0, -l0/6],
        [r0, 0, -l0/2],
    ])

    X_SiN = numpy.array([
        [Rx, 0, -(l0/2-lsin)],
        [rsin +tan*lsin, 0, -l0/2 +lsin],
        [rsin, 0, -l0/2],
    ])

    p_SAM = [Point(x, lcCenter) for x in X_SAM_ctr[:-1]]
    p_SAM.append(Point(X_SAM_ctr[-1],lcCenter/5.0))
    p_SiN = [Point(x, lcOuter) for x in X_SiN]
    p_Fluid_up = [Point(x, lcOuter) for x in X_Fluid_up]
    p_Fluid_low = [Point(x, lcOuter) for x in X_Fluid_low]
    p_Fluid_ctr = [Point(x, lcCenter) for x in X_Fluid_ctr]

    if sam is None:
        p_Au = []
    else:
        X_Au = numpy.array([
            [Rx, 0, -lsam +l0/2],
            [rsam -tan*(lsam-l0), 0, -lsam + l0/2],
            [rsam, 0, -l0/2],
        ])
        p_Au = [Point(X_Au[0], lcOuter)]
        p_Au.extend([Point(x, lcCenter) for x in X_Au[1:]])

    # Group all fluid points into list
    p_Fluid = p_Fluid_up + p_SAM
    if sam is not None:
        p_Fluid.append(p_Au[-1])