# -*- coding: utf-8 -*-


### Import Libraries
import os, tempfile
from pylab import *
import csv
from CSXCAD  import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
from pathlib import Path

import pandas as pd
import numpy as np

import skrf as rf

### Setup the simulation

Sim_Path = os.path.join(os.getcwd(), 'MSL20')
post_proc_only = False

unit = 1e-6 # specify everything in um
MSL_length = 8000
#MSL_length = 4000
MSL_width = 450
substrate_thickness = 203

substrate_width=8000
substrate_epr = 3.55
#stub_length = 12e3
f_max = 20e9

### Setup FDTD parameters & excitation function
FDTD = openEMS()
FDTD.SetGaussExcite( f_max/2, f_max/2 )
FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] )

### Setup Geometry & Mesh
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)

resolution = C0/(f_max*sqrt(substrate_epr))/unit/50 # resolution of lambda/50
third_mesh = array([2*resolution/3, -resolution/3])/4

## Do manual meshing
#mesh.AddLine('x', 0)
#mesh.AddLine('x',  MSL_width/2+third_mesh)
#mesh.AddLine('x', -MSL_width/2-third_mesh)
#mesh.SmoothMeshLines('x', resolution/2)

mesh.AddLine('x', 0)
#mesh.AddLine('x',  MSL_width/2+third_mesh)
#mesh.AddLine('x', -MSL_width/2-third_mesh)


mesh.AddLine('x', [-MSL_length/2, MSL_length/2])
mesh.SmoothMeshLines('x', resolution/2)

mesh.AddLine('y', 0)
mesh.AddLine('y',  MSL_width/2+third_mesh)
mesh.AddLine('y', -MSL_width/2-third_mesh)
mesh.SmoothMeshLines('y', resolution/4)

#mesh.AddLine('y', [-15*MSL_width, 15*MSL_width])

mesh.AddLine('y', [-substrate_width/2, substrate_width/2])
mesh.AddLine('y', (MSL_width/2))
mesh.SmoothMeshLines('y', resolution)

mesh.AddLine('z', linspace(0,substrate_thickness,5))
mesh.AddLine('z', 3000)
mesh.SmoothMeshLines('z', resolution)

## Add the substrate
substrate = CSX.AddMaterial( 'RO4350B', epsilon=substrate_epr)
start = [-MSL_length/2, -substrate_width/2, 0]
stop  = [+MSL_length/2, substrate_width/2, substrate_thickness]
substrate.AddBox(start, stop )

###Add the gnd plane
"""
gndp = CSX.AddMetal( 'PEC' )
startg = [-MSL_length, -15*MSL_width, 0]                                                                                            
stopg  = [ +MSL_length, +15*MSL_width, 0]
gndp.AddBox(startg, stopg, priority=10 )
"""
### add gnd as a conducting sheet

pec = CSX.AddConductingSheet('PEC', conductivity=59.6e6, thickness=35e-6)
startg = [-MSL_length/2, -substrate_width/2, 0]  
stopg  = [ +MSL_length/2, +substrate_width/2, 0]
pec.AddBox(startg, stopg, priority=10 )


## MSL Definition

topC = CSX.AddConductingSheet( 'copper' , conductivity=59.6e6, thickness=35e-6)
startp = [-MSL_length/2,  -MSL_width/2, substrate_thickness]                                                                                            
stopp  = [ MSL_length/2,  +MSL_width/2, substrate_thickness]
topC.AddBox(startp, stopp, priority=10 )

"""
Original 
pec = CSX.AddConductingSheet( 'PEC' , conductivity=59.6e6, thickness=35e-6)
startp = [-MSL_width/2,  MSL_width/2, substrate_thickness]                                                                                            
stopp  = [ MSL_width/2,  MSL_width/2, substrate_thickness]
pec.AddBox(startp, stopp, priority=10 )

"""
## MSL port setup
port = [None, None]

portstarta = [ -MSL_length/2, -MSL_width/2, substrate_thickness]
portstopa  = [ -MSL_length/4,  MSL_width/2, 0]
#port[0] = FDTD.AddMSLPort( 1,  pec, portstarta, portstopa, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=-MSL_length/6, priority=10)
port[0] = FDTD.AddMSLPort( 1,  topC, portstarta, portstopa, 'x', 'z', excite=-1,FeedShift=10*resolution, MeasPlaneShift=MSL_length/4, priority=10)
portstartb = [+MSL_length/2, -MSL_width/2, substrate_thickness]
portstopb  = [+MSL_length/4  ,  MSL_width/2, 0]
#port[1] = FDTD.AddMSLPort( 2, topC, portstartb, portstopb, 'x', 'z', MeasPlaneShift=MSL_length/4, priority=10 )
port[1] = FDTD.AddMSLPort( 2, topC, portstartb, portstopb, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=MSL_length/4, priority=10)


### Run the simulation
if 0:  # debugging only
    CSX_file = os.path.join(Sim_Path, 'MSL20.xml')
    if not os.path.exists(Sim_Path):
        os.mkdir(Sim_Path)
    CSX.Write2XML(CSX_file)
    from CSXCAD import AppCSXCAD_BIN
    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))


if not post_proc_only:
    FDTD.Run(Sim_Path, cleanup=True)
    
CSX_file = os.path.join(Sim_Path, 'MSL20.xml')
CSX.Write2XML(CSX_file)
from CSXCAD import AppCSXCAD_BIN
os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

freq = rf.Frequency(start=1, stop=20, npoints=1601, unit='ghz')
### Post-processing and plotting

#f = rf.Frequency(start=1, stop=10, npoints=101, unit='ghz')

#original
f = linspace( 1e9, f_max, 1601 )
for p in port:
    p.CalcPort( Sim_Path, f, ref_impedance = 50)

s11 = port[0].uf_ref / port[0].uf_inc
s21 = port[1].uf_ref / port[0].uf_inc

s22 = port[1].uf_ref / port[1].uf_inc
s12 = port[0].uf_ref / port[1].uf_inc


A = ((1+s11)*(1-s11) + s21*s21)/(2*s21);
B =  port[0].Z_ref*((1+s11)*(1-s11) - s21*s21)/(2*s21);

C = ((1-s11)*(1-s11) - s21*s21)/(2*s21) / port[1].Z_ref;
#C = ((1-s11)*(1-s22) - s21*s21)/(2*s21) / 50

D =  ((1+s11)*(1-s22) + s21*s12)/(2*s21);
#D=  ((1+s11)*(1-s22) + s21*s12)/(2*s21)/ 50

Zo_calc=sqrt(B/C);

s11_real=s11.real;
s11_imag=s11.imag;

s21_real=s21.real;
s21_imag=s21.imag;


s12_real=s12.real;
s12_imag=s12.imag;


s22_real=s22.real;
s22_imag=s22.imag;




# Create a plot
#plt.figure(figsize=(10, 6))  # Optional: set figure size


plt.plot(f, 20*log10(abs(s11)), 'k-', f, 20*log10(abs(s21)), 'r-', f,20*log10(abs(s22)), 'b-',f,20*log10(abs(s12)),'g-',linewidth = 8) # Example: line 

plt.legend(["S11", "S21", "S22", "S12"], loc="upper right")

# Add labels and title
plt.xlabel('Frequency',fontsize=20)
plt.ylabel('Sparam(dB)',fontsize=20)
plt.title('Sparams',fontsize=20)
plt.grid(True) # Optional: add a grid
plt.xlim(0, 20e9) # Set x-axis limit from 0 to 20
#plt.ylim(-45,100)
# Display the plot
plt.show()
        

