# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              
# 
# Path: sucemfem/Utilities/MeshGenerators.py
# def get_centred_cube(domain_size, max_edge_len, centred_element_coordinate=None):
#     """ Generate a cube mesh centred arount [0,0,0]
# 
#     optionally translates mesh slightly such that
#     centred_element_coordinate is at the centre of an element.
#     """
#     domain_subdivisions = N.array(N.ceil(N.sqrt(2)*domain_size/max_edge_len), N.uint)
#     mesh = dolfin.UnitCube(*domain_subdivisions)
#     # Transform mesh to correct dimensions
#     mesh.coordinates()[:] *= domain_size
#     mesh.coordinates()[:] -= domain_size/2
#     if centred_element_coordinate is not None:
#         ## Translate mesh slightly so that source coordinate lies at
#         ## centroid of an element
#         centred_element_point = dolfin.Point(*centred_element_coordinate)
# #        source_elnos = mesh.all_intersected_entities(centred_element_point)
#         source_elnos = mesh.intersected_cells(centred_element_point)
#         closest_elno = source_elnos[(N.argmin(
#             [centred_element_point.distance(dolfin.Cell(mesh, i).midpoint())
#              for i in source_elnos]))]
#         centre_pt = dolfin.Cell(mesh, closest_elno).midpoint()
#         centre_coord = N.array([centre_pt.x(), centre_pt.y(), centre_pt.z()])
#         # The mesh intersect operator caches certain data
#         # structures. We need to clear them if the mesh coordinates
#         # are changed after calling any mesh intersection methods.
#         mesh.intersection_operator().clear()
#         mesh.coordinates()[:] -= centre_coord
#     return mesh
# 
# Path: sucemfem/PostProcessing/surface_ntff.py
# class SurfaceNTFFForms(object):
# class NTFF(object):
#     def __init__(self, function_space):
#     def set_dofs(self, dofs):
#     def get_N_form(self):
#     def get_L_form(self):
#     def assemble_N(self):
#     def assemble_L(self):
#     def set_parms(self, theta, phi, k0):
#     def __init__(self, function_space):
#     def set_dofs(self, dofs):
#     def set_frequency(self, frequency):
#     def calc_pt(self, theta_deg, phi_deg):
#         N = self.get_N_form()
#         L = self.get_L_form()
# 
# Path: sucemfem/PostProcessing/variational_ntff.py
# class TransformTestingExpression(object):
# class NTFF(object):
#     def __init__(self):
#     def set_parms(self, rhat, ahat, k0):
#     def _set_expr_parms(self, expr, rhat, k0, constfac):
#     def get_expression(self):
#     def __init__(self, function_space, testing_space=None):
#     def set_k0(self, k0):
#     def set_frequency(self, frequency):
#     def set_dofs(self, dofs):
#     def calc_pt(self, theta_deg, phi_deg):
#     def calc_pt_E_H(self, theta_deg, phi_deg):
#     def calc_ff_func(self, rhat, ahat):


import pickle
import numpy as N
import sys
import sucemfem.Utilities.Optimization
import matplotlib as mpl
import matplotlib.pyplot as plt
from dolfin import *
from sucemfem.Consts import Z0, c0
from sucemfem.Utilities.MeshGenerators import get_centred_cube
from sucemfem.PostProcessing import surface_ntff, variational_ntff
V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", data['order'])

#------------------------------
# Calculate using the standard surface integral NTFF
surf_ntff = surface_ntff.NTFF(V)
surf_ntff.set_dofs(data['x'])
surf_ntff.set_frequency(data['freq'])
surf_E_ff = N.array([surf_ntff.calc_pt(th_deg, ph_deg)
                for th_deg, ph_deg in zip(theta_deg, phi_deg)])
surf_E_theta = surf_E_ff[:,0]
surf_E_phi = surf_E_ff[:,1]
#------------------------------
# Calculate using the variational integral NTFF
var_ntff = variational_ntff.NTFF(V)
var_ntff.set_k0(k0)
var_ntff.set_dofs(data['x'])
var_ntff.set_frequency(data['freq'])
var_E_ff = N.array([var_ntff.calc_pt(th_deg, ph_deg)
                for th_deg, ph_deg in zip(theta_deg, phi_deg)])
var_E_theta = var_E_ff[:,0]
var_E_phi = var_E_ff[:,1]
#------------------------------
# Analytical solution of the far-field of an infintesimal dipole
E_theta_an = Z0*N.exp(-1j*k0)*(1j*k0*1*lam/1000)/(4*N.pi)*N.sin(N.deg2rad(theta_deg))
#------------------------------
# Plot results
font = {'size'   : 18}

mpl.rc('font', **font)


# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              
# 
# Path: sucemfem/Utilities/MeshGenerators.py
# def get_centred_cube(domain_size, max_edge_len, centred_element_coordinate=None):
#     """ Generate a cube mesh centred arount [0,0,0]
# 
#     optionally translates mesh slightly such that
#     centred_element_coordinate is at the centre of an element.
#     """
#     domain_subdivisions = N.array(N.ceil(N.sqrt(2)*domain_size/max_edge_len), N.uint)
#     mesh = dolfin.UnitCube(*domain_subdivisions)
#     # Transform mesh to correct dimensions
#     mesh.coordinates()[:] *= domain_size
#     mesh.coordinates()[:] -= domain_size/2
#     if centred_element_coordinate is not None:
#         ## Translate mesh slightly so that source coordinate lies at
#         ## centroid of an element
#         centred_element_point = dolfin.Point(*centred_element_coordinate)
# #        source_elnos = mesh.all_intersected_entities(centred_element_point)
#         source_elnos = mesh.intersected_cells(centred_element_point)
#         closest_elno = source_elnos[(N.argmin(
#             [centred_element_point.distance(dolfin.Cell(mesh, i).midpoint())
#              for i in source_elnos]))]
#         centre_pt = dolfin.Cell(mesh, closest_elno).midpoint()
#         centre_coord = N.array([centre_pt.x(), centre_pt.y(), centre_pt.z()])
#         # The mesh intersect operator caches certain data
#         # structures. We need to clear them if the mesh coordinates
#         # are changed after calling any mesh intersection methods.
#         mesh.intersection_operator().clear()
#         mesh.coordinates()[:] -= centre_coord
#     return mesh
# 
# Path: sucemfem/PostProcessing/surface_ntff.py
# class SurfaceNTFFForms(object):
# class NTFF(object):
#     def __init__(self, function_space):
#     def set_dofs(self, dofs):
#     def get_N_form(self):
#     def get_L_form(self):
#     def assemble_N(self):
#     def assemble_L(self):
#     def set_parms(self, theta, phi, k0):
#     def __init__(self, function_space):
#     def set_dofs(self, dofs):
#     def set_frequency(self, frequency):
#     def calc_pt(self, theta_deg, phi_deg):
#         N = self.get_N_form()
#         L = self.get_L_form()
# 
# Path: sucemfem/PostProcessing/variational_ntff.py
# class TransformTestingExpression(object):
# class NTFF(object):
#     def __init__(self):
#     def set_parms(self, rhat, ahat, k0):
#     def _set_expr_parms(self, expr, rhat, k0, constfac):
#     def get_expression(self):
#     def __init__(self, function_space, testing_space=None):
#     def set_k0(self, k0):
#     def set_frequency(self, frequency):
#     def set_dofs(self, dofs):
#     def calc_pt(self, theta_deg, phi_deg):
#     def calc_pt_E_H(self, theta_deg, phi_deg):
#     def calc_ff_func(self, rhat, ahat):


# Path: examples/ntff/cube_ntff.py
import pickle
import numpy as N
import sys
import sucemfem.Utilities.Optimization
import matplotlib as mpl
import matplotlib.pyplot as plt
from dolfin import *
from sucemfem.Consts import Z0, c0
from sucemfem.Utilities.MeshGenerators import get_centred_cube
from sucemfem.PostProcessing import surface_ntff, variational_ntff

V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", data['order'])

#------------------------------
# Calculate using the standard surface integral NTFF
surf_ntff = surface_ntff.NTFF(V)
surf_ntff.set_dofs(data['x'])
surf_ntff.set_frequency(data['freq'])
surf_E_ff = N.array([surf_ntff.calc_pt(th_deg, ph_deg)
                for th_deg, ph_deg in zip(theta_deg, phi_deg)])
surf_E_theta = surf_E_ff[:,0]
surf_E_phi = surf_E_ff[:,1]
#------------------------------
# Calculate using the variational integral NTFF
var_ntff = variational_ntff.NTFF(V)
var_ntff.set_k0(k0)
var_ntff.set_dofs(data['x'])
var_ntff.set_frequency(data['freq'])
var_E_ff = N.array([var_ntff.calc_pt(th_deg, ph_deg)
                for th_deg, ph_deg in zip(theta_deg, phi_deg)])
var_E_theta = var_E_ff[:,0]
var_E_phi = var_E_ff[:,1]
#------------------------------
# Analytical solution of the far-field of an infintesimal dipole
E_theta_an = Z0*N.exp(-1j*k0)*(1j*k0*1*lam/1000)/(4*N.pi)*N.sin(N.deg2rad(theta_deg))
#------------------------------
# Plot results
font = {'size'   : 18}

mpl.rc('font', **font)


plt.figure()