# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/giov27/OpticalConductivityGraphene/blob/main/Store_OpticalConductivity.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + id="aqHzrgT7giPz"
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
import time
import os
# scipy.__version__
# from numba import guvectorize, complex64
# from scipy import constants

# print(constants.Boltzmann)
# For local runtime
# jupyter notebook --NotebookApp.allow_origin='https://colab.research.google.com' --port=8888 --NotebookApp.port_retries=0 --no-browser
# pd.set_option("max_rows", None)
# pd.set_option("max_colwidth", None)

# + colab={"base_uri": "https://localhost:8080/"} id="lOz8bFKDtz2T" outputId="0f856552-e0c6-4913-a6ca-e67c5cb30087"
os.getcwd()

# + id="L6cyHrcB313p"
# from google.colab import drive
# drive.mount('/content/drive', force_remount=True)

# + colab={"base_uri": "https://localhost:8080/"} id="Mtvb7tGAgodC" outputId="c5fe64fb-2248-49b6-fe29-571b9bc56785"
# constanta

na = 2048 # Number of atom in both direction
nm = na-1
nx = na # Number of atom in x-direction
ny = na # Number of atom in y-direction
nl = na//2 # Number of repetation for 2x2 block matrix
lx = 1 # Length of wavefunction in x-direction
ly = 1 # Length of wavefunction in y-direction
dx = lx/nx # Interval x
dy = ly/ny # Interval y

nrx = nx//2
nry = ny//2

pot = 0 # On-site potensial 

nt = 4096 # timestep
nta = nt//2 # Rounded half of timestep
deltaT = 0.02 # interval time
ta = (nt/2)*deltaT # Maximum value of Time Propagation
Tm = ta*(np.pi/(2*np.pi))
dTm = Tm/nta

a = 1 # hopping parameter
mu = 0 # Chemical Potential

T = 300 # temperature
# Kb = 0.03564 #1 # boltzmann constant
Kb = 0.695 #1 # boltzmann constant
# T_au = T /(3.158*10**(5))
beta = 1 / (Kb * T)

sum0 = 0

print(dTm)
print(Tm)

# + id="Y5C86kgfqh_P"
# Array Declaration
H_even = np.zeros((nx,ny), dtype=float)
H_odd = np.zeros((nx,ny), dtype=float)
fermiDirac = np.zeros((nx,ny), dtype=complex)
fermiDirac_t1 = np.zeros((nx,ny), dtype=complex)
# fermiDirac_t2 = np.zeros((nx,ny), dtype=complex)

# + id="3YmhRV8dkh7f"
## Declaration of array Hamiltonan
def H(H, num): # num = 0 for even, num = 1 for odd
  if num == 0:
    d = 0
    for i in range(nrx):
      for j in range(nry):
        if j == i:
          sx = 2*i+d
          sy = 2*j+d
          H[sx][sy+1] = -a
          H[sx+1][sy] = -a

  if num == 1:
    H[0][nm] = -a
    H[nm][0] = -a
    d = 1
    for i in range(nrx-1):
      for j in range(nry-1):
        if j == i:
          sx = 2*i+d
          sy = 2*j+d
          H[sx][sy+1] = -a
          H[sx+1][sy] = -a
  return H

## Normalize Hamiltonan
def normalize_H(H, num, sum0):
  sum0n = 0
  for i in range(nx):
    for j in range(ny):
      sum0n += H[i][j]*H[i][j]
      # print(sum0)

  if num == 0:
    sum0 += np.sqrt(4*sum0n)
  elif num == 1:
    sum0 += np.sqrt(2*sum0n)
  return sum0


# + id="omFfdsten92w"
def fermi(ff, c, psi, mat = fermiDirac, mat_t1 = fermiDirac_t1): 
  # ff = 1 for 1-FermiDirac, c related to normal(0), conjugate(1)
  # MFD = 0 # Matrix Fermi Dirac
  TEn = Kb*T # Bottom part of fermi dirac distribution
  potN = pot/sum0 # Normalize potential
  bN = sum0/TEn # Normalize beta (BetaN = Beta*HamiltonanNorm )
  aN = a/sum0 # Normalize hopping parameter

  ch = np.cosh(a/TEn) # Cos hiperbolic part
  sh = np.sinh(bN*aN) # Sin hiperbolic part

  ## Filling the inverse matrix of fermi dirac: Diagonal Part
  for i in range(nx):
    for j in range(ny):
      if j == i:
        mat[i][j] = ch * np.exp((bN*potN) - (mu/TEn))

  ## Filling the inverse matrix of fermi dirac: 2x2 Block diagonal Part with even index
  d = 0
  for i in range(nrx):
    for j in range(nry):
      if i == j :
        sx = 2*i+d
        sy = 2*j+d
        mat[sx][sy+1] = -sh * np.exp((bN*potN) - (mu/TEn))
        mat[sx+1][sy] = -sh * np.exp((bN*potN) - (mu/TEn))
   
  ## Filling the inverse matrix of fermi dirac: 2x2 Block diagonal with odd index
  d = 0
  i = 0
  for j in range(ny):
    sx = 2*i+d
    sy = j
    mat_t1[sx][sy] = ch*mat[sx][sy] - sh*mat[nm][sy]
    mat_t1[nm][sy] = -sh*mat[sx][sy] + ch*mat[nm][sy]

  d = 1
  for i in range(nrx-1):
    for j in range(ny):
      sx = 2*i+d
      sy = j
      mat_t1[sx][sy] = ch*mat[sx][sy] - sh*mat[sx+1][sy]
      mat_t1[sx+1][sy] = -sh*mat[sx][sy] + ch*mat[sx+1][sy]

  ## Adding identity matrix to mat_t
  mat_t2 = np.identity(nx) + mat_t1

  ## Inverse Matrix
  invMat = np.linalg.inv(mat_t2)

  if ff == 0:
    if c == 0:
      invMat_t = invMat
    elif c == 1:
      invMat_t = np.conjugate(invMat)

  elif ff == 1:
    invMat_t = np.identity(nx) - invMat

  ## Multiplication (1-F(H-mu)^-1) with psi0
  psi1_t = invMat_t @ psi

  faktorNormalisasi = 1/(np.linalg.norm(psi1_t)*np.sqrt(dx*dy))
  psi1_t *= faktorNormalisasi

  return psi1_t


# + id="A48l1nEC_YuG"
def currentDens(psi1):
  psi = np.zeros((nx,ny), dtype = complex)
  psi_a = np.zeros((nx,ny), dtype = complex)
  c0 = 0 + 0j
  s0 = 0 - 1j*a
  # s0 = 0 - 0.5j*a
  d=1
  for i in range(nl-1):
    for j in range(ny):
      sx = 2*i+d
      sy = j
      sx1 = sx+1
      psi[sx][sy] = c0*psi1[sx][sy] + s0*psi1[sx1][sy]
      psi[sx1][sy] = s0*psi1[sx][sy] + c0*psi1[sx1][sy]

  i = 0
  for j in range(ny):
    sx = 2*i
    sy = j
    sx1 = sx+1
    psi[sx][sy] =  c0*psi1[sx][sy] + s0*psi1[nm][sy]
    psi[nm][sy] =  s0*psi1[sx][sy] + c0*psi1[nm][sy]
  
  d=1
  for i in range(nrx):
    for j in range(nry):
      sx = 2*i+d
      sy = 2*j+d
      sx = (sx+1) % nx
      psi_a[sx][sy] = c0*psi1[sx][sy] + s0*psi1[sx1][sy]
      psi_a[sx1][sy] = s0*psi1[sx][sy] + c0*psi1[sx1][sy]
  
  d=0
  for i in range(nrx):
    for j in range(nry):
      sx = 2*i+d
      sy = 2*j+d
      sx1 = sx+1
      psi_a[sx][sy] = c0*psi1[sx][sy] + s0*psi1[sx1][sy]
      psi_a[sx1][sy] = s0*psi1[sx][sy] + c0*psi1[sx1][sy]

  psi0 = psi + psi_a
  
  faktorNormalisasi = 1/(np.linalg.norm(psi0)*np.sqrt(dx*dy))
  psi0 *= faktorNormalisasi
  return psi0


# + id="7IavThMlzq3D"
def kineticPropagation(d, t, psi1, deltaT):
  er = np.cos(-deltaT*a)
  ei = -1.0j*(np.sin(-1.0*deltaT*a))
  # er = np.cos(-1.0*dTm*a)
  # ei = -1.0j*(np.sin(-1.0*dTm*a))

  if(d == 0):
    nry = int(ny/2)
    for x in range(nx):
      for y in range(nry):
        sx = x
        sy = int((2*y) + t)
        sy1 = int((sy+1)%ny)
        psi1[sx][sy] = er*psi1[sx][sy] - ei*psi1[sx][sy1]
        psi1[sx][sy1] = -ei*psi1[sx][sy] + er*psi1[sx][sy1]

  elif(d == 1):
    nrx = int(nx/2)
    nry = int(ny/2)
    for x in range(nrx):
      for y in range(nry):
        sx = (2*x) + t
        sy = (2*y) + t
        sx1 = (sx+1)%nx
        psi1[sx][sy] = er*psi1[sx][sy] - ei*psi1[sx1][sy]
        psi1[sx1][sy] = -ei*psi1[sx][sy] + er*psi1[sx1][sy]
  return psi1


# + id="Q55TDspCQHH2"
## One Cycle of kinetic propagation
def oneCycle(psi, deltaT):
  for d in range(2):
    for t in range(2):
      kineticPropagation(d, t, psi, deltaT)

  faktorNormalisasi = 1/(np.linalg.norm(psi)*np.sqrt(dx*dy))
  psi *= faktorNormalisasi
  return psi


# + id="hNKTHvIM6FL3"
## Declare psi1 
def psi1_evol():
  psi1_t = np.zeros((nta, nx, ny), dtype=complex)
  psi1_t[0] = fermiDirac1_t
  # print('{} = {}'.format(0, psi1_t[0]))
  for k in range(1,nta):
    psi1_t1 = oneCycle(psi1_t[k-1], deltaT)
    psi1_t[k] = psi1_t1
    # print('{} = {}'.format(k, psi1_t1))
  return psi1_t


# + id="jVrx9G63wce7"
## Declare psi2
def psi2_evol():
  psi2_t = np.zeros((nt,nx,ny), dtype=complex)
  psi2_t[0] = fermiDirac2_t
  # print('{} = {}'.format(0, psi2_t[0]))
  # y= input()
  for k in range(1,nt):
    psi2_t1 = oneCycle(psi2_t[k-1], deltaT)
    psi2_t[k] = psi2_t1
    # print('{} = {}'.format(k, psi2_t1))
    # y= input()
  return psi2_t


# + id="BTrwiLLDcddF"
def omega():
  # Tm = ta*(1/(2*np.pi)) #new periode is shifted by factor 0.159(or 1/(2*pi)) hence the freq window following Nyquist Theorm
  # dTm = (Tm/nta) #delta Tau maks
  print('Tau Maks: {}, Delta Tau: {}'.format(Tm, dTm))
  tau = 0
  tau_arr = []
  for k in range(nt):
    if (tau <= Tm):
      tau_arr.append(tau)
      tau = tau + dTm

  #Constanta Opcon
  Wm = 1/(2*dTm) # Maximum frequency, with Nyquist Theorem applied
  dWm = 1/(2*Tm) # Frequency interval, with Nyquist Theorem applied
  print('W Maks: {}, Delta W: {}'.format(Wm, dWm))
  W = 0.001
  W_arr = np.linspace(0, Wm, nta)
  W_arr2 = []
  for k in range(nt):
    if (W <= Wm):
      W_arr2.append(W)
      W = W + dWm
    
  return tau_arr, W_arr,W_arr2,Tm, Wm


# + colab={"base_uri": "https://localhost:8080/"} id="sAzxMlyCFEHC" outputId="3d98e81e-481f-4f47-d8b1-6f1ff6943ea2"
tau_arr, W_arr,W_arr2, tau_max, W_max = omega()

# + colab={"base_uri": "https://localhost:8080/"} id="GvjSYvZiFazZ" outputId="e90cb1c6-94de-48c3-d443-33b038cabf0f"
print(W_arr)

# + colab={"base_uri": "https://localhost:8080/"} id="DMQbtJrFGTJp" outputId="9a730592-d8cc-4b11-cdcc-8c970a2671e4"
print(W_arr2)


# + id="eIyG3pM26xrC"
## Correlation of Optical Conductivity
def opCon_corr():
  cfr_arr = []
  # psi1_t = np.array(psi1_evol())
  psi1_t = psi1_evol()
  # psi2_t = np.conjugate(np.array(psi2_evol()))
  psi2_t = psi2_evol()
  opCon = []
  cfr_arr = []

  for k in range(nt):
    J = np.array(currentDens(psi1_t[k]))

    opCon_t = psi2_t[k] @ J

    opCon.append(opCon_t)
    # opCon.append(2*(opCon_t).imag)
  opCon_t2 = np.array(opCon)

  # multiplication opCon[0] with opCon[k]
  for k in range(nt):
    cfr = 0
    for i in range(nx):
      for j in range(ny):
        cfr += opCon_t2[0][i][j] * opCon_t2[k][i][j]
    # cfr_arr.append((cfr*dx*dy).imag)
    cfr_arr.append(cfr*dx*dy)

  return np.array(cfr_arr)


# + id="MaQPUFqZqrzy"
psi0 = np.random.uniform(-1,1,(nx,ny)) + 1j*np.random.uniform(-1,1,(nx,ny))
faktorNormalisasi = 1/(np.linalg.norm(psi0)*np.sqrt(dx*dy))
psi0 *= faktorNormalisasi

psi1 = psi0
psi2 = psi0.conjugate()

# pd.DataFrame(psi0)

# + colab={"base_uri": "https://localhost:8080/"} id="jxP9uLJKc_6i" outputId="54e9835a-0c2a-4087-9139-18394d04bf83"
H(H_odd, 1)
# pd.DataFrame(H(H_odd, 1))

# + colab={"base_uri": "https://localhost:8080/"} id="DHUAd-zzc6EK" outputId="f52eb2b4-efea-4d1c-bed2-145c7e88b695"
H(H_even, 0)
# pd.DataFrame(H(H_even, 0))

# + colab={"base_uri": "https://localhost:8080/"} id="nrDp8PDgdIk5" outputId="f77b15ca-214c-47ba-ca7f-57dfb5501a2a"
sum0 = normalize_H(H_even, 0, sum0)
sum0 = normalize_H(H_odd, 1, sum0)
print(sum0)


# + id="h_s_JTXEXbhL"
## Correlation of Opctical Conductivity with combine psi
def opCon_corr2(psi1):
  opCon = []
  opCon_t = np.zeros((nt,nx,ny), dtype=complex)
  cfr_arr = np.zeros((nt), dtype=complex)

  psi1_J = currentDens(psi1)
  fermiDirac1_t = fermi(1, 0, psi1_J)
  psi1_J2 = currentDens(fermiDirac1_t)
  fermiDirac2_t = fermi(0, 1, psi1_J2)

  for k in range(nt):
    if k == 0:
      opCon_t[0] = fermiDirac2_t
    else :
      kinProp1 = oneCycle(opCon_t[k-1], 2*deltaT)
      # kinProp2 = oneCycle(kinProp1, deltaT)
      opCon_t[k] = kinProp1

    # multiplication psi1[0] with opCon[k]
    cfr = 0
    for i in range(nx):
      for j in range(ny):
        cfr += psi2[i][j] * opCon_t[k][i][j]
    # cfr_arr.append((cfr*dx*dy).imag)
    cfr_arr[k]= cfr*dx*dy
    # cfr_arr[k]= 2*cfr*dx*dy

  return cfr_arr


# + id="1JdJcx6ei8mE"
# try to implement the saving memory with odd even login opcon, and save the correlation
# def corr_func(k, opCon_t):
def corr_func(opCon_t):
  cfr = 0
  for i in range(nx):
    for j in range(ny):
      cfr += psi2[i][j] * opCon_t[i][j]
      # cfr += psi2[i][j] * opCon_t[k][i][j]
  # cfr_arr.append((cfr*dx*dy).imag)
  # cfr_arr[k]= cfr*dx*dy
  return cfr*dx*dy

def opCon_corr3(psi1):
  # opCon = []
  opCon_t = np.zeros((3,nx,ny), dtype=complex)
  cfr_arr = np.zeros((nta), dtype=complex)

  psi1_J = currentDens(psi1)
  fermiDirac1_t = fermi(1, 0, psi1_J)
  psi1_J2 = currentDens(fermiDirac1_t)
  # fermiDirac2_t = fermi(0, 1, psi1_J2)
  fermiDirac2_t = fermi(0, 0, psi1_J2)

  for k in range(nta):
    if k == 0:
      opCon_t[0] = fermiDirac2_t
      cfr_arr[k] = corr_func(opCon_t[0])      

    elif k == 1:
      kinProp1 = oneCycle(opCon_t[0], 2*dTm)
      opCon_t[1] = kinProp1
      cfr_arr[k] = corr_func(opCon_t[1])  
    
    elif k%2 == 0:
      kinProp1 = oneCycle(opCon_t[1], 2*dTm)
      opCon_t[2] = kinProp1
      cfr_arr[k] = corr_func(opCon_t[2])

    else :
      kinProp1 = oneCycle(opCon_t[2], 2*dTm)
      # kinProp2 = oneCycle(kinProp1, deltaT)
      opCon_t[1] = kinProp1
      cfr_arr[k] = corr_func(opCon_t)  

    # cfr_arr[k]= cfr*dx*dy
    # cfr_arr[k]= 2*cfr*dx*dy

  return cfr_arr


# + id="R2bostHyYZrw"
# try to implement the saving memory with odd even login opcon, and save the correlation
def opCon_corr3(psi1):
  # opCon = []
  opCon_t = np.zeros((3,nx,ny), dtype=complex)
  cfr_arr = np.zeros((nta), dtype=complex)
  tau = 0

  psi1_J = currentDens(psi1)
  fermiDirac1_t = fermi(1, 0, psi1_J)
  psi1_J2 = currentDens(fermiDirac1_t)
  # fermiDirac2_t = fermi(0, 1, psi1_J2)
  fermiDirac2_t = fermi(0, 0, psi1_J2)

  for k in range(nta):
    if tau<Tm or k<=1000:
      if k == 0:
        opCon_t[0] = fermiDirac2_t
        cfr = corr_func(opCon_t[0])
        # cfr_arr[k] = 1*np.exp(-1*0.18*tau)*cfr
        cfr_arr[k] = 2*np.exp(-1*0.01*tau)*cfr


      elif k == 1:
        tic = time.time()
        kinProp1 = oneCycle(opCon_t[0], 2*deltaT)
        opCon_t[1] = kinProp1
        cfr = corr_func(opCon_t[1])
        cfr_arr[k] = 2*np.exp(-1*0.01*tau)*cfr
        toc = time.time()
        running_time = toc - tic
        print(running_time)
      
      elif k%2 == 0:
        kinProp1 = oneCycle(opCon_t[1], 2*deltaT)
        opCon_t[2] = kinProp1
        cfr = corr_func(opCon_t[2])
        cfr_arr[k] = 2*np.exp(-1*0.01*tau)*cfr
        

      else :
        kinProp1 = oneCycle(opCon_t[2], 2*deltaT)
        # kinProp2 = oneCycle(kinProp1, deltaT)
        opCon_t[1] = kinProp1
        cfr = corr_func(opCon_t[1])
        cfr_arr[k] = 2*np.exp(-1*0.01*tau)*cfr  

      # cfr_arr[k]= cfr*dx*dy
      tau = dTm*(k)
      print('iterasi:', k)


    # cfr_arr[k]= 2*cfr*dx*dy

  return cfr_arr

# + id="bs8eBrdHZ2-i" colab={"base_uri": "https://localhost:8080/", "height": 442} outputId="f1ca9a67-9ffa-4e77-f32a-5cce71a45daa"
#at the start of cell, define <tic> variable
tic = time.time()

opCon_corr2_psi3 = opCon_corr3(psi0)

#at the end of cell, define <toc> variable
toc = time.time()

running_time = toc - tic
print(running_time)

# + id="Ov97Dk8-me9E"
# tic = time.time()

# opCon_corr2_psi2 = opCon_corr2(psi0)

#at the end of cell, define <toc> variable
# toc = time.time()

# running_time = toc - tic
# print(running_time)

# + id="hc0zSCN5nMrv"
print(opCon_corr2_psi3[3])
# print(opCon_corr2_psi2[3])

# + id="5pe-fOCnmtfY"
plt.plot(np.arange(0, nta), opCon_corr2_psi3, label='3 part')
# plt.plot(np.arange(0, nt), opCon_corr2_psi2, label='2 part')
# plt.xlim(-1, 1000)
plt.legend()
plt.show()

# + id="ZDpSdLhmQd3z"
header_save = 'Optical Conductivity {} x {} \n nt={} ; deltaT={}; Tm={}; dTm={} \n Kb={} ; T={}; time={} s'.format(na,na,nt,deltaT,Tm, dTm, Kb,T,running_time)

# + id="FBd1lSI-yJXG"
# pd.DataFrame(opCon_corr2_psi0)

# + id="CN_e0e6ax6sk"
print(opCon_corr2_psi3.shape)

# + id="KG76_zoPzAuw"
# np.savetxt('/content/drive/MyDrive/Data/opCon1.csv', opCon_corr2_psi0.T,fmt='%.12e%+.12ej', delimiter=",", header=header_save,footer=header_save, comments='~~~ ')

#true
# np.savetxt('/content/drive/MyDrive/Data/Data3/opCon{}_{}.csv'.format(na,Kb), opCon_corr2_psi3.T,fmt='%.12e%+.12ej', delimiter=",",  header=header_save, footer=header_save)
np.savetxt('opCon{}_{}.csv'.format(na,Kb), opCon_corr2_psi3.T,fmt='%.12e%+.12ej', delimiter=",",  header=header_save, footer=header_save)

# + id="AZTwP_Lw2UaM"
# # !cat '/content/drive/MyDrive/Data/Data2/opCon32_0.695.csv'
# # !cat '/opCon512_0.695.csv'

# + id="EXxVuzdh29TO"
# data = np.loadtxt('/content/drive/MyDrive/Data/Data2/opCon{}_{}.csv'.format(na,Kb), delimiter=',' , dtype=complex)
data = np.loadtxt('opCon{}_{}.csv'.format(na,Kb), delimiter=',' , dtype=complex)
pd.DataFrame(data)

# + id="v_YyN-bwr5vV"
# psi1_J = currentDens(psi1)
# fermiDirac1_t = fermi(1, 0, psi1_J)
# fermiDirac2_t = fermi(0, 0, psi2)

# pd.DataFrame(fermiDirac1_t)

# + id="TTET0Get2EAe"
# print(psi1_evol())

# + id="xh-Rc_E82YeS"
# print(psi2_evol())

# + id="MU56SlyqRFmm"
# opCon_t2 = opCon_corr()

# + id="nANfBCvWPCN6"
# x = [elem.real for elem in opCon_t2]
# y = [elem.imag for elem in opCon_t2]
x = [elem.real for elem in opCon_corr2_psi3]
y = [elem.imag for elem in opCon_corr2_psi3]
# x = [elem.real for elem in data]
# y = [elem.imag for elem in data]

# + id="BP41cBF0R5Na"
plt.plot(np.arange(0, nta), x, label='real part')
plt.plot(np.arange(0, nta), y, label='imag part')
# plt.xlim(-1, 1000)
plt.legend()
plt.show()

# + id="rSBTks1o0wnz"
tau_arr, W_arr,W_arr2, tau_max, W_max = omega()

# + id="f0gg432iktPt"
# opCon_t3 = opCon_calc()
# opCon_t4 = opCon_calc2(tau_arr)
# opCon_t4 = opCon_t2
# print(opCon_t4)

# + id="k80MIeDz7Kh3"
print(len(tau_arr))
print(len(opCon_corr2_psi3))
tau_arr1 = np.linspace(-tau_max, tau_max, nt)
print(tau_arr1)

# + id="-OaopQdT12Pt"
X = np.zeros(len(opCon_corr2_psi3))
fungsiKorelasi = np.hstack((X,opCon_corr2_psi3))
print(len(fungsiKorelasi))
# X = np.zeros(len(opCon_corr2_psi0))
# fungsiKorelasi = np.hstack((X,opCon_corr2_psi0))

plt.plot(tau_arr1, fungsiKorelasi.real, label='real part')
plt.plot(tau_arr1, fungsiKorelasi.imag, label='imag part')
# plt.xlim(0, 6)
# plt.legend()

# + id="Ov2VIo8o9Ucd"
# plt.plot(np.arange(-len(opCon_t2), len(opCon_t2)), fungsiKorelasi.real, label='real part')
# plt.plot(tau_arr1, fungsiKorelasi.real, label='real part')
# plt.xlim(-1, 1)
# plt.legend()

# + id="pQbWDpD69pnZ"
# plt.plot(tau_arr1, fungsiKorelasi.imag, label='Imag part')
# plt.xlim(-2, 2)
# plt.legend()

# + id="YrGRgoFA0nu3"
##Opcon FFT numpy

def opConFFT(func):
  hann = np.hanning(len(func))
  # wp = func.real*hann
  # wp = func.imag*hann
  wp = func*hann

  y = np.fft.fft(wp)/(2*np.pi)
  # y = np.fft.rfft(wp)/(2*np.pi)
  # y1 = np.abs(y)
  y1 = y
  # y1 = y.imag
  
  # freq = np.fft.fftfreq(len(func), 0.95*deltaT)
  freq = np.fft.fftfreq(len(func), 5.9*dTm)
  # freq = np.fft.rfftfreq(len(func),0.095)
  # print(freq)
  freq1 = np.fft.fftshift(freq*(2*np.pi)/1)
  # print(freq1)

  result = np.fft.fftshift(y1)
  print(len(wp))
  # result = np.abs(y1)
  return freq1, result


# + id="HtQXnapY9JSS"
opConFFT_t = opConFFT(fungsiKorelasi)
# opConFFT_t = opConFFT(opCon_t2)
plt.title('Optical Condictivity\n[2 Dimension]')
plt.plot(opConFFT_t[0], opConFFT_t[1])
plt.xlim(0,10)
plt.grid()
plt.show()

# + id="Wm6drtI-6DVz"
plt.plot(opConFFT_t[0], np.absolute(opConFFT_t[1].imag), label='Imag part')
# plt.plot(opConFFT_t[0], np.absolute(opConFFT_t[1].real), label='Real part')
# plt.plot(opConFFT_t[0], np.absolute(opConFFT_t[1]), label='Both part')
plt.xlim(0,15)
plt.legend()
plt.grid()
plt.show()


# + id="E_d-RUnH3U8Z"
def constApply(func):
  xMid = len(func[0])//2
  xLast = len(func[0])
  Wm = 1/(2*dTm) # Maximum frequency, with Nyquist Theorem applied
  dWm = 1/(2*ta)
  opCon_x = func[0][xMid:xLast]
  print(len(opCon_x))
  omega = np.linspace(0,6,xMid)
  opCon_valueSlice = func[1][xMid:xLast]
  # const = -beta*omega
  # const2 = (np.exp(const)-1)/omega
  # print(const2)
  # print(const)
  opCon_Const = []
  opCon = []
  for k in opCon_x:
    if k == 0 :
      opCon_Const.append(0)
    else :
      const1 = -beta*k
      opCon_Const.append(-(np.exp(const1)-1)/k)
  # for k in range(xMid, xLast,1):
  #     print(k)
  #     tau = dTm*(k+2)
  #   if k == 0 :
  #     opCon_Const.append(0)
  #   else :
  #     const1 = -beta*k
  #     # print(np.exp(const1))
  #     opCon_Const.append((np.exp(const1)-1)/k)
  const = np.array(opCon_Const)
  for k in range(xMid):
    # opCon.append(((const[k] * func[1][k].imag)/1))
    # opCon.append(((const[k]*200 * opCon_valueSlice[k])/1))
    opCon.append(((const[k]* 134.33*opCon_valueSlice[k].imag)/1))
    # opCon.append(((opCon_Const[k] * func[1][k])/0.04).real)
  # opCon_t6 = np.fft.fftshift(opCon)
  print(opCon_Const)
  return opCon_x, opCon


# + id="JANAlQto49pG"
newValue = constApply(opConFFT_t)

# + id="3Jq8o7YI5KXc"
plt.plot(newValue[0],np.abs(newValue[1]))
plt.xlim(0,10)


# + id="uE2IttHgG-Ij"
## opcon FFT scipy

def opConFFT_scipy(func):
  g = np.hanning(len(func))
  # wp = np.abs(func*g)
  wp = func.imag*g

  y = scipy.fft.fft(wp)/(2*np.pi)
  y1 = y
  result = scipy.fft.fftshift(y1)

  freq = scipy.fft.fftfreq(len(func), deltaT)
  # print(freq)
  # freq1 = scipy.fft.fftshift(freq*((2*np.pi)/3.05))
  freq1 = scipy.fft.fftshift(freq*((2*np.pi)))

  return freq1, result


# + id="tUaithwAJd9w"
# opConFFT_sci = opConFFT_scipy(fungsiKorelasi)
# max_value = np.max(opConFFT_sci[1].imag)
# print(4/max_value)
# plt.title('Optical Condictivity\n[2 Dimension]')
# plt.plot(opConFFT_sci[0], np.abs(opConFFT_sci[1]))
# plt.plot(opConFFT_sci[0], np.abs(opConFFT_sci[1]))
# plt.xlim(-16,16)
# plt.grid()
# plt.show()

# + id="mzxpIJV9QjRD"
# plt.plot(opConFFT_sci[0], opConFFT_sci[1].imag, label='Imag part')
# plt.plot(opConFFT_sci[0], opConFFT_sci[1].real, label='Real part')
# plt.xlim(-7,7)
# plt.legend()
# plt.grid()
# plt.show()

# + id="Aa3VE7eSQsEt"
# plt.title('Optical Condictivity\n[2 Dimension]')
# plt.plot(opConFFT_sci[0], opConFFT_sci[1].imag)
# plt.xlim(0,8)
# plt.grid()
# plt.show()

# + id="6N-xJoSVpcVn"
# result = np.where(opConFFT_t[0] == 0)
# print(result)

# + id="73u1cWhT5R7k"
def constApply():
  xVal = len(opConFFT_t[0])
  opCon_Const = []
  opCon = []
  for k in opConFFT_t[0]:
    if k == 0 :
      opCon_Const.append(0)
    else :
      const1 = -beta*k
      # print(np.exp(const1))
      opCon_Const.append((np.exp(const1)-1)/k)
  const = np.array(opCon_Const)
  for k in range(xVal):
    opCon.append(((const[k] * opConFFT_t[1][k].imag)/1))
    # opCon.append(((opCon_Const[k] * opConFFT_t[1][k])/0.04).real)
  # opCon_t6 = np.fft.fftshift(opCon)
  return opCon

# + id="OmmFXsnp2ooT"
# plt.plot(opConFFT_t[0], np.abs(opCon_t6))
# plt.xlim(-0.7,10)
# plt.ylim(0, 2.5)
# plt.grid()
# plt.show()

# + id="vI3p1VnKumMq"
# plt.plot(opConFFT_t[0], opCon_t6)
# plt.xlim(-1,10)
# plt.ylim(-0.0002, 0.0002)
# plt.grid()
# plt.show()

StarcoderdataJupyterScriptsDedupFiltered