setup_helper_functions.py

import pandas as pd
import numpy as np
import scipy as sci


# # # # # set_nf_subsystems(simP, Fock_matrix) # # # # #

# One method of saving all cases of fock_trunc info for each qubit or cavity mode considered in the dressed eigenstate of the system
# you can write another method for your purpose

# Outputs: all different cases of fock_trunc info in the qubits and cavity modes
# in each rows of the matrix

# inputs an array of fock_truc to uniformly apply to all the qubits in the
# system
# inputs an array of fock_trunc to uniformly apply to all the cavity modes
# in the system

# inputs:
# simP
# Fock_matrix

# outputs: the updated simP
# simP['nf_qubits']
# simP['nf_cavity']
# simP['num_cases']


def set_nf_subsystems(simP, Fock_matrix):
    if np.shape(Fock_matrix)[1] != simP['num_mode']:
        raise ValueError('put column size of the Fock_matrix as equal to the total number of qubits and '
                         'cavity modes in the dressed eigenstate of the system')

    simP['nf_qubits'] = Fock_matrix[:,
                               0:simP['num_t']]  # 02/04/2024 changed name from simP.nf_qubits_cases
    simP['nf_cavity'] = Fock_matrix[:,
                               simP['num_t']:simP['num_mode']]  # 02/04/2024 from simP.nf_cavity_cases
    return simP

# end set_nf_subsystems()

# # # # # rect_freq_cavity_perturbation(simP, mode_string, probe_len) # # # # #
def rect_freq_cavity_perturbation(simP, mode_string, probe_len):

    mode_chr = mode_string[:2]
    m = float(mode_string[2])
    n = float(mode_string[3])
    p = float(mode_string[4])

    num_rm = sum([simP['reg'][i]['reg_Type'] == 0 for i in range(len(simP['reg']))])  # for kk modes
    num_cm = sum([simP['reg'][i]['reg_Type'] == 1 for i in range(len(simP['reg']))])  # number of probes in the system
    f_unpert = np.array([])
    del_f = np.zeros((num_rm, num_cm))
    f_reson = np.zeros(num_rm)
    omega_Moon = np.zeros(num_rm)
    for aa in range(num_rm):

        # rectangular cavity's dimension by region
        a = simP['reg'][aa]['x']
        b = simP['reg'][aa]['y']
        d = simP['reg'][aa]['z']

        # cavity volume
        V_cavity = a * b * d

        # f_unpert analytical solution
        M = np.divide(m, a)  # a should be in m
        N = np.divide(n, b)  # b should be in m
        P = np.divide(p, d)  # d should be in m

        f_unpert = np.append(f_unpert, 0.5 / np.sqrt(simP['mu_0'] * simP['reg'][aa]['eps']) *
                             np.sqrt(np.power(M, 2) + np.power(N, 2) + np.power(P, 2))) # f_unpert = 0.5*(3*10^8)/(mu_r*eps_r_rm)*sqrt(M.^2+N.^2+P.^2);

        # assuming k is determined by the geometry of the cavity
        k = np.sqrt(np.power(m * np.pi / a, 2) + np.power(n * np.pi/ b, 2) + np.power(p * np.pi / d, 2))

        for bb in range(num_cm):

            rc = simP['reg'][num_rm+bb]['inner_radius']  # rc = simP.reg(2).inner_radius;

            # probes
            V_probe = np.pi * rc**2 * probe_len
            # probe location
            x0 = simP['reg'][num_rm+bb]['x0']
            y0 = simP['reg'][num_rm+bb]['y0'](probe_len)
            z0 = simP['reg'][num_rm+bb]['z0']

            if mode_chr == 'TE':
                if n == 0:
                    mode_factor = 4
                elif m == 0:
                    mode_factor = 4
                else:
                    mode_factor = 8

                # Field expressions
                Hx_c = (m * np.pi / a) * (p * np.pi / d)
                Hx = Hx_c * np.sin(m * np.pi / a * x0) * np.cos(n * np.pi / np.multiply(b, y0)) * np.cos(p * np.pi / d * z0)

                Hy_c = (n * np.pi / b) * (p * np.pi / d)
                Hy = Hy_c * np.cos(m * np.pi / a * x0) * np.sin(n * np.pi / np.multiply(b, y0)) * np.cos(p * np.pi / d * z0)

                ktmn = np.sqrt((m * np.pi / a)**2 + (n * np.pi / b)**2)
                Hz = ktmn**2 * np.cos(m * np.pi / a * x0) * np.cos(n * np.pi / np.multiply(b, y0)) *\
                     np.sin(p * np.pi / d * z0)

                Ex_c = (n * np.pi / b)
                Ex = Ex_c * np.cos(m * np.pi / a * x0) * np.sin(n * np.pi / np.multiply(b, y0)) * np.sin(p * np.pi / d * z0)

                Ey_c = (m * np.pi / a)
                Ey = Ey_c * np.sin(m * np.pi / a * x0) * np.cos(n * np.pi / np.multiply(b, y0)) * np.sin(p * np.pi / d * z0)

                # perturbation ratio eqaution expression
                numer_c = (np.multiply(Hx, 2) + np.multiply(Hy, 2) + np.multiply(Hz, 2)) * simP['mu_0'] -\
                          (np.multiply(Ex, 2) + np.multiply(Ey, 2)) *\
                          (k * simP['reg'][aa]['eta'])**2 * simP['reg'][aa]['eps']
                numer = np.multiply(numer_c, V_probe)

                denom_c = (Hx_c**2 + Hy_c**2 + ktmn**4) * simP['mu_0'] + (Ex_c**2+Ey_c**2) *\
                          (k * simP['reg'][aa]['eta'])**2 * simP['reg'][aa]['eps']
                denom = denom_c * V_cavity / mode_factor



            if mode_chr == 'TM':

                # mode factor
                if p == 0:
                    mode_factor = 4
                else:
                    mode_factor = 8


                # Field expressions
                Hx_c = (n * np.pi / b)
                Hx = Hx_c * np.sin(m * np.pi / a * x0) * np.cos(n * np.pi / np.multiply(b, y0)) *\
                     np.cos(p * np.pi / d * z0)

                Hy_c = (m * np.pi / a)
                Hy = Hy_c * np.cos(m * np.pi / a * x0) * np.sin(n * np.pi / np.multiply(b, y0)) *\
                     np.cos(p * np.pi / d * z0)

                Ex_c = (m * np.pi / a) * (p * np.pi / d)
                Ex = Ex_c * np.cos(m * np.pi / a * x0) * np.sin(n * np.pi / np.multiply(b, y0)) *\
                     np.sin(p * np.pi / d * z0)

                Ey_c =  (n * np.pi / b) * (p * np.pi / d)
                Ey = Ey_c * np.sin(m * np.pi / a * x0) * np.cos(n * np.pi / np.multiply(b, y0)) *\
                     np.sin(p * np.pi / d * z0)

                ktmn = np.sqrt((m * np.pi / a)**2 + (n * np.pi / b)**2)
                Ez = ktmn**2 * np.sin(m * np.pi / a * x0) * np.sin(n * np.pi / np.multiply(b, y0)) *\
                     np.cos(p * np.pi / d * z0)

                # perturbation ratio eqaution expression
                numer_c = (np.multiply(Hx, 2) + np.multiply(Hy, 2)) * (k / simP['reg'][aa]['eta'])**2 *\
                          simP['mu_0'] - (np.multiply(Ex, 2) + np.multiply(Ey, 2) + np.multiply(Ez, 2)) *\
                          simP['reg'][aa]['eps']
                numer = np.multiply(numer_c, V_probe)

                denom_c = (Hx_c**2 + Hy_c**2) * (k / simP['reg'][aa]['eta'])**2 * simP['mu_0'] +\
                          (Ex_c**2 + Ey_c**2 + ktmn**4) * simP['reg'][aa]['eps']
                denom = denom_c * V_cavity / mode_factor

            if probe_len == 0:
                ratio = 0  # At zero perturbation.
            else:
                # common equation
                ratio = numer / denom


            del_f[aa][bb] = np.multiply(ratio, f_unpert[aa])
        f_reson[aa] = f_unpert[aa] + sum(del_f[aa])  # bug caught 11/21/2023 f_reson(aa) = f_unpert(aa) + sum(del_f(aa,bb),2);
        omega_Moon[aa] = 2 * np.pi * f_reson[aa]
    return omega_Moon, f_reson
# End rect_freq_cavity_perturbation()

# # # # # cal_cavity_eigfreq(simP) # # # # #

# calls the subfunction "rect_freq_cavity_perturbation" to calculate the ports probes in cavity's resonant
# frequencies using the cavity perturbation theory
# and outputs the real part of the resonant frequencies
# both linear and angular parts

# input
# simP
# simP['cavity_mode_string']
# simP['probe_len']
#
# output simP
# simP['f_cavity_arr']
# simP['omega_r']

def cal_cavity_eigfreq(simP):
    simP['f_cavity_arr'] = np.zeros(simP['num_r'])
    omega_Moon1 = np.zeros(simP['num_r'])
    f_cavity_res1 = np.zeros(simP['num_r'])
    for ir in range(simP['num_r']):
        # outputs: cavity resonant frequency in linear and angular
        temp = rect_freq_cavity_perturbation(simP, simP['cavity_mode_string'][ir], simP['probe_len'])  # 7.5525 GHz
        omega_Moon1[ir] = temp[0]
        f_cavity_res1[ir] = temp[1]

    # assuming one cavity region
    for ir in range(simP['num_r']):
        # outputs
        simP['f_cavity_arr'][ir] = np.real(f_cavity_res1[ir]) # (1): means (the region 1)% only consider the real part
        simP['omega_r'] = 2 * np.pi * simP['f_cavity_arr']  # rad/s
    # more cavity region then, additional code for that needs to be written
    return simP
# End cal_cavity_eigfreq()

# # # Capacitance_of_dipole_ant() # # #

# calculates the analytical capacitance of the dipole antenna
# using the Balanis textbook formula
# assumes 90 degrees orientation

# input: simP
# simP['dL']
# simP['f_cavity_arr']
# simP['c_freespace']
# simP['a_radius']

# output: simP
# simP['Cap_Bal']

def Capacitance_of_dipole_ant(simP):

    dL_Bal = simP['dL']  # in m
    Cap_Bal = np.zeros(1)
    for ir in range(1):  #:simP['num_r'] # calculates capacitance in respect to the fundamental cavity mode % comment added 02/10/2024
        # parameter values for analytically computing the capacitance of dipole

        f_Bal = simP['f_cavity_arr'][ir]  # Hz
        omega_Bal = 2 * np.pi * f_Bal
        k_Bal = omega_Bal / simP['c_freespace']

        # expresions for the computated constants
        # Balanis Antenna Theory 4th edition
        # small dipole antenna input reactance & capacitance analytical expression
        # Balanis eqn. 8-59
        Xm = -120 * np.divide((np.log(dL_Bal / (2 * simP['a_radius'])) - 1), (np.tan(k_Bal * dL_Bal / 2)))  # in ohms

        # outputs simP.Cap_Bal
        Cap_Bal[ir] = -(np.divide(1, (omega_Bal * Xm)))
    # output: the capacitance of the dipole antenna at different resonant
    # frequencies
    simP['Cap_Bal'] = Cap_Bal  # in F
    return simP
# End Capacitance_of_dipole_ant()

# # # rect_cavity_field() # # #

def rect_cavity_field(*argv): # (simP,polarP, aa,bb, mode) %x,y,z,eta,string_mode,m,n,p,a,b,d) % x,y,z in m
    # if (nargin == 8)
    #     [simP,polarP, aa,bb, mode] = deal(varargin{4:end});
    Field = dict()
    if (len(argv) == 6):
        simP = argv[3]
        cavity = argv[4]
        mode = argv[5]
        #         x = argv[0];
        #         y = argv[1];
        #         z = argv[2];
    elif (len(argv) == 4):
        simP = argv[0]
        polarP = argv[1]
        cavity = argv[2]
        mode = argv[3]

    # cavity dimension
    a = cavity['x']
    b = cavity['y']
    d = cavity['z']

    string_mode = mode[:2]
    m = float(mode[2])
    n = float(mode[3])
    p = float(mode[4])

    if (len(argv) == 6):
        x = argv[0]
        y = argv[1]
        z = argv[2]
    elif (len(argv) == 4):  # (nargin == 5)
        # rename the variable
        x = polarP['xq3']  # xqq
        y = b
        z = polarP['yq3']  #zq

    # Analytical field expressions
    # assuming k is determined by the geometry of the cavity
    k = np.sqrt((m * np.pi / a)**2 + (n * np.pi / b)**2 + (p * np.pi / d)**2)
    # k_mnp = sqrt(  (m*pi/a)^2 + (n*pi/b)^2  + (p*pi/d)^2); % added for the condensed form     06/07/2023
    ktmn = np.sqrt((m * np.pi / a)**2 + (n * np.pi / b)**2)
    eps_r = cavity['eps'] / simP['eps_0']

    if string_mode[:2] == 'TE':

        #mode factor
        if n == 0:
            mode_factor = 4
        elif m == 0:
            mode_factor = 4
        else:
            mode_factor = 8


        # for TE mode
        # compute the normalization constant
        Hmnp = 1
        # sqrt(N_u)
        N_u = np.sqrt(eps_r * (Hmnp * (k * cavity['eta']))**2 / ktmn**2 * (a * b * d / mode_factor)) # used to have negative sign inside the sqrt term 06/07/2023
        N_v = np.sqrt(Hmnp**2 / ktmn**2 * k**2 * (a * b * d / mode_factor))  # condensed form 06/07/2023

        # field equations
        Ex_c = Hmnp * ((-1) * (k * cavity['eta']) / ktmn**2) * (n * np.pi / b)  # replace the 1i term with the (-1)
        Field['Ex'] = (1 / N_u) * Ex_c * np.multiply(np.multiply(np.cos(m * np.pi / a * x),
                                                                 np.sin(n * np.pi / np.multiply(b, y))),
                                                     np.sin(p * np.pi / d * z))

        Ey_c = -Hmnp * ((-1) * (k * cavity['eta']) / ktmn**2) * (m * np.pi / a)  # replace the 1i term with the (-1)
        Field['Ey'] = (1 / N_u) * Ey_c * np.multiply(np.multiply(np.sin(m * np.pi / a * x),
                                                     np.cos(n * np.pi/ np.multiply(b, y))),
                                                     np.sin(p * np.pi / d * z))

        Field['Ez'] = np.zeros(np.size(x))  # Ez = zeros(x_row,x_col);% TE mode

        Hx_c = -Hmnp * (1 / ktmn**2) * (m * np.pi / a) * (p * np.pi / d)
        Field['Hx'] = (1 / N_v) * Hx_c * np.multiply(np.multiply(np.sin(m * np.pi / a * x),
                                                     np.cos(n * np.pi / np.multiply(b, y))), np.cos(p * np.pi / d * z))

        Hy_c = -Hmnp * (1 / ktmn**2) * (n * np.pi / b) * (p * np.pi / d)
        Field['Hy'] = (1 / N_v) * Hy_c * np.multiply(np.multiply(np.cos(m * np.pi / a * x),
                                                     np.sin(n * np.pi / np.multiply(b, y))), np.cos(p * np.pi / d * z))

        Field['Hz'] = (1 / N_v) * Hmnp * np.multiply(np.multiply(np.cos(m * np.pi / a * x),
                                                  np.cos(n * np.pi / np.multiply(b, y))),
                                                  np.sin(p * np.pi / d * z))
    if string_mode[:2] == 'TM':
        # mode factor
        if p == 0:
            mode_factor = 4
        else:
            mode_factor = 8

        # for TM mode
        # compute the normalization constant
        Emnp = 1

        N_u = np.sqrt(eps_r * Emnp**2 / ktmn**2 * k**2 * (a * b * d / mode_factor))
        N_v = np.sqrt((Emnp * (k / cavity['eta']))**2 / ktmn**2 * (a * b * d / mode_factor)) # used to have negative sign 06/07/2023

        # field equations
        Ex_c = -Emnp / ktmn**2 * (m * np.pi / a) * (p * np.pi / d)
        Field['Ex'] = (1 / N_u) * Ex_c * np.multiply(np.multiply(np.cos(m * np.pi / a * x),
                                                                 np.sin(n * np.pi / np.multiply(b, y))),
                                                     np.sin(p * np.pi / d * z))

        Ey_c = -Emnp / ktmn**2 * (n * np.pi / b) * (p * np.pi / d)
        Field['Ey'] = (1 / N_u) * Ey_c * np.multiply(np.multiply(np.sin(m * np.pi / a * x),
                                                  np.cos(n * np.pi / b * y)), np.sin(p * np.pi / d * z))

        Field['Ez'] = (1 / N_u) * Emnp * np.multiply(np.multiply(np.sin(m * np.pi / a * x),
                                                              np.sin(n * np.pi / b * y)),
                                                  np.cos(p * np.pi / d * z))

        Hx_c = Emnp * ((-1) * (k / cavity['eta']) / ktmn**2) * (n * np.pi / b) # replace the 1i term with the (-1)
        Field['Hx'] = (1 / N_v) * Hx_c * np.multiply(np.multiply(np.sin(m * np.pi / a * x),
                                                                 np.cos(n * np.pi / b * y)),
                                                     np.cos(p * np.pi / d * z))

        Hy_c = -Emnp * ((-1) * (k / cavity['eta']) / ktmn**2) * (m * np.pi / a)  # replace the 1i term with the (-1)
        Field['Hy'] = (1 / N_v) * Hy_c * np.multiply(np.multiply(np.cos(m * np.pi / a * x),
                                                                 np.sin(n * np.pi / b * y)),
                                                     np.cos(p * np.pi / d * z))

        Field['Hz'] = np.zeros(np.shape(x))

    return Field
# End rect_cavity_field()

# # # calculate_normVter(simP) # # #

# access each cavity subregion and analytically calculate the coupling strength
# between the cavity eigenmode and each qubit
# using the classical EM theory (vector effective length of a dipole antenna)

# the induced voltage across the unloaded dipole antenna (the open circuit voltage) by the cavity field
# is calculated by computing the incident normalized E-field at the location of the dipole inside the cavity
# using the subfunction "rect_cavity_field" and the vector effective length of the antenna.

# The normalized terminal voltage across the transmon (capacitively loaded dipole antenna in our formalism)
# is calculated by taking the real part of the open circuit voltage
# and using the voltage division to compute the voltage at the load (the Josephson Junction)

# input: simP
# simP['dL']
# simP['theta']
# simP['reg'][aa]: cavity subregion structure
# simP['cavity_mode_string']: kinds of modes that can exist in cavity
# simP['dipole'][it]['x0'], simP['dipole'][it]['y0'], simP['dipole'][it]['z0']: dipole location
# simP['Cap_Bal'], simP['Cap_load']
#
# output
# simP['normVter']

def calculate_normVter(simP):
    num_rm = sum([simP['reg'][i]['reg_Type'] == 0 for i in range(len(simP['reg']))])
    # vector effective length of a dipole oriented to an angle of theta from
    # the incident plane wave
    vector_effective_length = simP['dL'] / 2 * np.sin(simP['theta'])
    Vterminal = np.zeros((simP['num_t'], simP['num_r']))

    for aa in range(num_rm):  # for each cavity region
        cavity = simP['reg'][aa]

        for ir in range(simP['num_r']):  # for the considered number of cavity modes
            mode = simP['cavity_mode_string'][ir]

            for it in range(simP['num_t']):  # for all the qubits inside

                # dipole location
                x = simP['dipole'][it]['x0']
                y = simP['dipole'][it]['y0']
                z = simP['dipole'][it]['z0']

                Cavity_field = rect_cavity_field(x, y, z, simP, cavity, mode)
                Einc_of_cur_mode = Cavity_field['Ey']


                # Vterminal compuation
                init_Voc_analy = np.real(Einc_of_cur_mode * vector_effective_length)
                Vterminal[it][ir] = np.divide(np.multiply(init_Voc_analy, simP['Cap_Bal']), simP['Cap_Bal'] + simP['Cap_load'])

    simP['normVter'] = Vterminal

    return simP
# End calculate_normVter()