try.py

'''Convert the same into a function
freq step value changed to 0.01 instead of 0.001
'''
import scipy.stats as stats
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np
from numpy import pi, sqrt as isqrt, cos, sin
import random
import time
import pandas as pd

start = time.time()
print("Execution started at: ", dt.datetime.now())


fig, axs = plt.subplots(20, 4)
plt.figure(dpi=1200)


def calc_MA(da, kzA, pA):
    MA = []
    for i in range(10):
        MA.append(np.array(
            [[cos(kzA * da[i]), (1j / pA) * sin(kzA * da[i])], [1j * pA * sin(kzA * da[i]), cos(kzA * da[i])]]))
    return MA


def calc_MB(db, kzB, pB):
    return np.array([[cos(kzB * db), (1j / pB) * sin(kzB * db)], [1j * pB * sin(kzB * db), cos(kzB * db)]])  # VO2


def calc_MC(dc, kzC, pC):
    MC_def = []
    for i in range(10):
        MC_def.append(np.array(
            [[cos(kzC * dc[i]), (1j / pC) * sin(kzC * dc[i])], [1j * pC * sin(kzC * dc[i]), cos(kzC * dc[i])]]))
    return MC_def


c1 = 0
c2 = 0

finalR, finalAs, finalTs, finalAp, finalTp = [[], [], [], [], []]


def plotUsingSigma(x, y, sigma=200, variance=0.0):
    global c1, c2

    st = dt.datetime.now()
    print("Starting iteration: ", int(x / 10), " ", y, "time: ", dt.datetime.now())

    # plt.rcParams["figure.figsize"] = [7.00, 3.50]
    # plt.rcParams["figure.autolayout"] = True

    # Constants
    c = 3e8
    e0 = 8.85e-12
    u0 = (4 * pi) * (10 ** -7)
    e = 1.6e-19
    h = 6.63e-34
    hbar = h / (2 * pi)

    # Constants 2
    n0 = 1  # air
    nt = 1
    ut = 1
    et = nt ** 2  # t:substrate
    na = 2.25
    ua = 1
    ea = na ** 2  # SiO2
    nc = 3.3
    uc = 1
    ec = nc ** 2  # Si
    ub = 1
    lambda0 = 70e-6
    sigma_VO2 = sigma  # S/m Variable 200 or 200000

    # *****************************VO2
    eps_inf = 12
    gamma = 5.75e13  # rad/s
    wp = 1.4e15 * isqrt(sigma_VO2 / 3e5)
    theta = 0 * pi / 180  # Variable
    TETM = 1
    jj = 0  # set to 0 as python indices start from 0
    n = 0
    # Special Variables

    # da = lambda0 / (4 * na)  # SiO2
    # dc = lambda0 / (4 * nc)  # Si
    # db = lambda0 / (4 * sqrt(12))  # VO2

    # variance = 1.0e-6
    mean = 6.25e-6
    data_points = 10

    # MA, MC, R, Ts, As, Tp, Ap, FRE = [[], [], [], [], [], [], [], []]

    for i in range(data_points):

        # da = list(np.random.normal(mean, variance, data_points))  # find truncated
        # db = lambda0 / (4 * isqrt(12))  # VO2
        # dc = list(np.random.normal(mean, variance, data_points))
        da = stats.truncnorm.rvs(((1e-6 - mean) / variance), (11.25e-6 - mean) / variance, loc=mean, scale=variance,
                                 size=10)
        db = lambda0 / (4 * isqrt(12))  # VO2
        dc = stats.truncnorm.rvs(((1e-6 - mean) / variance), (11.25e-6 - mean) / variance, loc=mean, scale=variance,
                                 size=10)
        MA, MC, R, Ts, As, Tp, Ap, FRE = [[], [], [], [], [], [], [], []]
        jj = 0
        for fre in np.arange(start=0.1, stop=20.001, step=0.01) * 1e12:  # orig: step 0.001
            w = 2 * pi * fre
            epsVO2 = eps_inf - wp ** 2 / (w ** 2 + 1j * gamma * w)  # used j instead of i
            ###
            kz0 = (w / c) * n0 * cos(theta)
            kzA = (w / c) * na * isqrt(1 - (sin(theta) ** 2 / (na ** 2)))
            kzB = (w / c) * (isqrt(epsVO2)) * (
                isqrt(1 - (sin(theta) ** 2 / epsVO2)))  # isqrt for sqrt of complex number
            kzC = (w / c) * nc * isqrt(1 - (sin(theta) ** 2 / (nc ** 2)))
            kzt = (w / c) * nt * isqrt(1 - (sin(theta) ** 2 / (et * ut)))

            if TETM == 1:
                y0 = u0
                ya = ua * u0
                yb = ub * u0
                yc = uc * u0
                yt = ut * u0  # TE
            else:
                y0 = -e0
                ya = -ea * e0
                yb = -epsVO2 * e0
                yc = -ec * e0
                yt = -et * e0  # TM

            p0 = kz0 / (w * y0)
            pA = kzA / (w * ya)
            pB = kzB / (w * yb)
            pC = kzC / (w * yc)
            pt = kzt / (w * yt)

            m1 = 0.5 * np.array([[1, -1 / p0], [1, 1 / p0]])  # 2x2 Matrix multiplied with scalar
            mt = np.array([[1], [-pt]])
            ###
            MA = calc_MA(da, kzA, pA)
            MB = calc_MB(db, kzB, pB)
            MC = calc_MC(dc, kzC, pC)
            ###
            X = np.linalg.multi_dot(
                [MA[0], MC[0], MA[1], MC[1], MA[2], MC[2], MA[3], MC[3], MA[4], MC[4], MB, MC[5], MA[5], MC[6],
                 MA[6],
                 MC[7], MA[7], MC[8], MA[8], MC[9], MA[9]])
            ###
            MM = np.linalg.multi_dot([m1, X, mt])
            r = MM[1][0] / MM[0][0]
            t = 1 / (MM[0][0])
            ar = float(abs(r) ** 2)

            R.append(ar)
            ###
            Ts.append((kzt / kz0) * abs(t) ** 2)
            As.append(1 - R[jj] - Ts[jj])
            ###
            Tp.append((kzt / (et * kz0)) * abs(t) ** 2)
            Ap.append(1 - R[jj] - Tp[jj])
            ###
            FRE.append(fre)
            jj = jj + 1

        axs[i + x, y].plot(FRE, Ts, color='green', linestyle='dotted')  # Transmittance
        axs[i + x, y].plot(FRE, R, color='black', linestyle='dashed')  # Reflectance
        axs[i + x, y].plot(FRE, As, color='blue', linestyle='solid')  # Absorption
        axs[i + x, y].set_xticks([])
        axs[i + x, y].set_yticks([])
        finalR.append(R)
        finalTp.append(Tp)
        finalAp.append(Ap)
        finalTs.append(Ts)
        finalAs.append(As)

        c1 += 1
    print("Took : ", dt.datetime.now() - st)


variance_list = [0.5e-6, 1.0e-6, 1.5e-6, 2.0e-6]
sigma_values = [200, 200000]
a = 0
b = 0

for s in sigma_values:
    for v in variance_list:
        if sigma_values.index(s) == 0:
            x = 0
        else:
            x = 10
        plotUsingSigma(x, variance_list.index(v), sigma=s, variance=v)
plt.autoscale(enable=True)
plt.figure(dpi=600)

plt.show()

end = time.time()
print("Execution ended at: ", dt.datetime.now())
print("Execution Time: ", (end - start), "")