# -*- coding: utf-8 -*-
# input data

import numpy as np
#import matplotlib.pyplot as plt

def read(fn):
    # initialize
    Title = ''
    Nfreq = 0
    Freq = []
    Iplanewave = 0
    Planewave = [0, 0, 0, 0]
    Iradiusall = 2
    Radiusall = 0.2
    Iground = 0
    Z0 = 50
    Plot3dgeom = 0

    # read file
    with open(fn, 'rt', encoding='utf-8') as fp:
        header = fp.readline()
        lines = fp.readlines()

    # header
    d = header.strip().split()
    if (len(d) < 3) or (d[0] != 'OpenMOM'):
        print('*** Not OpenMOM data.')
        return
    #version = (100 * int(d[1])) + int(d[2])
    
    # number of geometries
    ngeom = 0
    for tline in lines:
        #print(tline)
        if tline.strip().startswith('geometry = '):
            d = tline.strip().split()
            ngeom += 6 if (d[2] == '2') and ((d[3] == '21') or (d[3] == '22') or (d[3] == '23')) else 1
    Ngeom = ngeom
    #print(Ngeom)
    
    # alloc geometry data
    G_gtype   = np.zeros( ngeom,        int)
    G_cosys   = np.zeros( ngeom,        int)
    G_pos     = np.zeros((ngeom, 4, 3), float)
    G_div     = np.zeros((ngeom, 2),    int)
    G_ifeed   = np.zeros( ngeom,        int)
    G_feed    = np.zeros((ngeom, 2),    float)
    G_iload   = np.zeros( ngeom,        int)
    G_load    = np.zeros( ngeom,        float)
    G_iradius = np.zeros( ngeom,        int)
    G_radius  = np.zeros( ngeom,        float)
    G_offset  = np.zeros((ngeom, 3),    float)
    
    # body
    ng = -1
    for tline in lines:
        # token
        d = tline.strip().split()  #;print(d)
        
        # check
        if d[0].startswith('#'):
            continue
        elif (d[0] == 'end'):
            break
        elif (len(d) < 3):
            continue
        elif (d[1] != "="):
            continue

        # parameter array (string)
        p = d[2:]  #;print(p)
        
        # input
        key = d[0]
        if   key == 'title':
            #Parm['title'] = tline[len('title = ') :].strip()
            Title = tline[len('title = ') :].strip()
        elif key == 'geometry':
            if len(p) > 5:
                ngeom, gtype, cosys, pos, div = geometry_data(p)
                #print(ngeom, gtype, cosys, pos, div)
                for igeom in range(ngeom):
                    ng += 1
                    G_gtype[ng] = gtype[igeom]
                    G_cosys[ng] = cosys[igeom]
                    G_pos[ng]   = pos[igeom]
                    G_div[ng]   = div[igeom]
        elif key == 'feed':
            if len(p) > 1:
                G_ifeed[ng] = 1
                G_feed[ng] = [float(p[0]), float(p[1])]  # V,deg
        elif key == 'load':
            if len(p) > 1:
                G_iload[ng] = int(p[0])  # = 1/2/3
                G_load[ng] = float(p[1])   # = R/L/C
        elif key == 'radius':
            G_iradius[ng] = 1
            G_radius[ng] = float(p[0])
        elif key == 'offset':
            if len(p) > 2:
                for ig in range(ng - ngeom + 1, ng + 1):
                    #Geom[ng - ngeom + 1 : ng]["offset"] = [p[0], p[1], p[2]]
                    G_offset[ig] = [float(p[0]), float(p[1]), float(p[2])]
        elif key == 'name':
            pass
            #Geom[ng]["name"] = tline[len('name = ') :]
        elif key == 'planewave':
            if len(p) > 2:
                Iplanewave = int(p[2])
                Planewave[0:2] = [float(p[0]), float(p[1])]
                if (int(p[2]) == 5) and (len(p) > 4):
                    Planewave[2] = float(p[3])
                    Planewave[3] = float(p[4])
        elif key == 'ground':
            Iground = int(p[0])
        elif key == 'frequency':
            if len(p) > 2:
                Freq = np.linspace(float(p[0]), float(p[1]), int(p[2]) + 1)
                Nfreq = len(Freq)
        elif key == 'radiusall':
            if len(p) > 1:
                Iradiusall = int(p[0])
                Radiusall = float(p[1])
        elif key == 'z0':
            Z0 = float(p[0])
        elif key == 'plot3dgeom':
            Plot3dgeom = int(p[0])
        #elif key == 'geom3dnode':
        #    if len(p) > 2:
        #        geom3dnode = [float(p[0]), int(p[1]), int(p[2])]
        #elif key == 'marker':
        #    marker = float(p[0])
        #elif key == 'window2d':
        #    if len(p) > 3:
        #        w2d = [int(p[0]), int(p[1]), int(p[2]), int(p[3])]
        #elif key == 'window3d':
        #    if len(p) > 5:
        #        w3d = [int(p[0]), int(p[1]), float(p[2]), int(p[3]), float(p[4]), float(p[5])]

    # font (3D)
    #plt.rcParams['font.family'] = 'sans-serif' if w3d[3] == 0 else ('serif' if w3d[3] == 1 else 'monospace')
    #plt.rcParams['font.size'] = w3d[2]
    #plt.rcParams['font.family'] = 'sans-serif' if w3d[3] == 0 else ('serif' if w3d[3] == 1 else 'monospace')
    #plt.rcParams['font.size'] = w3d[2]

    return Ngeom, G_gtype, G_cosys, G_pos, G_div, G_ifeed, G_feed, G_iload, G_load, G_iradius, G_radius, G_offset, \
        Nfreq, Freq, Title, Iplanewave, Planewave, Iground, Iradiusall, Radiusall, Z0, Plot3dgeom

# set geometry data
def geometry_data(p):
    g = int(p[0])
    s = int(p[1])
    
    ngeom = 6 if (g == 2) and ((s == 21) or (s == 22) or (s == 23)) else 1
    #ngeom = 1
    gtype = np.zeros( ngeom,        int)
    cosys = np.zeros( ngeom,        int)
    pos   = np.zeros((ngeom, 4, 3), float)
    div   = np.zeros((ngeom, 2),    int)
    #print(p)
    
    if (g == 1):
        #print(g, s)
		# wire unit
        if (s == 1) or (s == 2) or (s == 3):
            # point (1-2)
            gtype[0] = g
            cosys[0] = s
            pos[0, 0:4, 0] = [float(p[2]), float(p[3]), 0, 0] # X1, X2
            pos[0, 0:4, 1] = [float(p[4]), float(p[5]), 0, 0] # Y1, Y2
            pos[0, 0:4, 2] = [float(p[6]), float(p[7]), 0, 0] # Z1, Z2
            div[0, 0:2] = [p[8], 0]
        elif (s == 11):
            # X-line
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 0] = [float(p[2]), float(p[3]), 0, 0]  # X1, X2
            pos[0, 0:4, 1] = [float(p[4]), float(p[4]), 0, 0]  # Y
            pos[0, 0:4, 2] = [float(p[5]), float(p[5]), 0, 0]  # Z
            div[0, 0:2] = [p[6], 0]
        elif (s == 12):
            # Y-line
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 1] = [float(p[2]), float(p[3]), 0, 0]  # Y1, Y2
            pos[0, 0:4, 2] = [float(p[4]), float(p[4]), 0, 0]  # Z
            pos[0, 0:4, 0] = [float(p[5]), float(p[5]), 0, 0]  # X
            div[0, 0:2] = [p[6], 0]
        elif (s == 13):
            # Z-line
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 2] = [float(p[2]), float(p[3]), 0, 0]  # Z1, Z2
            pos[0, 0:4, 0] = [float(p[4]), float(p[4]), 0, 0]  # X
            pos[0, 0:4, 1] = [float(p[5]), float(p[5]), 0, 0]  # Y
            div[0, 0:2] = [p[6], 0]
    elif (g == 2):
        # plane unit
        if (s == 1) or (s == 2) or (s == 3):
            gtype[0] = g
            cosys[0] = s
            pos[0, 0:4, 0] = p[2:6]    # X1,X2,X3,X4
            pos[0, 0:4, 1] = p[6:10]   # Y1,Y2,Y3,Y4
            pos[0, 0:4, 2] = p[10:14]  # Z1,Z2,Z3,Z4
            div[0, 0:2] = p[14:16]     # N12,N14
            #print(pos, div)
        elif (s == 11):
            # X-rect
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 0] = p[2]
            pos[0, 0:4, 1] = [p[3], p[4], p[4], p[3]]
            pos[0, 0:4, 2] = [p[5], p[5], p[6], p[6]]
            div[0, 0:2] = p[7:9]
        elif (s == 12):
            # Y-rect
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 1] = p[2]
            pos[0, 0:4, 2] = [p[3], p[4], p[4], p[3]]
            pos[0, 0:4, 0] = [p[5], p[5], p[6], p[6]]
            div[0, 0:2] = p[7:9]
        elif (s == 13):
            # Z-rect
            gtype[0] = g
            cosys[0] = 1
            pos[0, 0:4, 2] = p[2]
            pos[0, 0:4, 0] = [p[3], p[4], p[4], p[3]]
            pos[0, 0:4, 1] = [p[5], p[5], p[6], p[6]]
            div[0, 0:2] = p[7:9]
        elif (s == 21) or (s == 22) or (s == 23):
            # box
            #ngeom = 6
            x1 = p[2]; x2 = p[3]
            y1 = p[4]; y2 = p[5]
            z1 = p[6]; z2 = p[7]
            nx = p[8]; ny = p[9]; nz = p[10]
            for m in range(ngeom):
                gtype[m] = g
                cosys[m] = s % 10
            # 1/2 : -X/+X
            for m in range(0, 2):
                pos[m, 0:4, 0] = x1 if (m == 0) else x2
                pos[m, 0:4, 1] = [y1, y2, y2, y1]
                pos[m, 0:4, 2] = [z1, z1, z2, z2]
                div[m, 0:2] = [ny, nz]
            # 3/4 : -Y/+Y
            for m in range(2, 4):
                pos[m, 0:4, 1] = y1 if (m == 2) else y2
                pos[m, 0:4, 2] = [z1, z2, z2, z1]
                pos[m, 0:4, 0] = [x1, x1, x2, x2]
                div[m, 0:2] = [nz, nx]
            # 5/6 : -Z/+Z
            for m in range(4, 6):
                pos[m, 0:4, 2] = z1 if (m == 4) else z2
                pos[m, 0:4, 0] = [x1, x2, x2, x1]
                pos[m, 0:4, 1] = [y1, y1, y2, y2]
                div[m, 0:2] = [nx, ny]
    
    return ngeom, gtype, cosys, pos, div
