# -*- coding: utf-8 -*- """ Created on Wed Oct 4 13:08:55 2017 @author: pwb504 """ import numpy as np import matplotlib.pyplot as plt import csv from scipy.optimize import curve_fit plt.close('all') ####### Constants and parameters ########################################### Energy_Array = [] min_freq = 50E6 max_freq = 3E9 # Only applies correction to signal between the freq limits Coil_Area = 9E-6 # Prodyn B-24 B-dot probe effective area Ddot_Area = 1.54E-5 # Equivalent area of D-dot FD-5C probe const = (3E8) / (2 * 1.26E-6) # c / (2 * mu0) ###### Cable Length ####### cable_length = 30 # length of cable in metres between probe and oscilloscope ###### Cable attenuation data for RG223 coax ######## # A list of frequencies and corresponding attenuation factors define the attenuation curve freq_MHz = [3.70, 150, 450, 824, 960, 1500, 2000, 2300, 5800] freq_interp = [i * 1E6 for i in freq_MHz] # convert from MHz to Hz so we can apply the correction properly att_dB_100ft = [-3.70, -8.10, -14.0, -19.0, -21.0, -26.7, -31.5, -34.1, -59.5] ##### Function to fit to cable attenuation data ##### def attenuation(f, A, B, C): return A * np.sqrt(f) + B + C*f #### Fitting and verification #### popt, pcov = curve_fit(attenuation, freq_interp, att_dB_100ft) # Fitting ######### A second cable attenuation function #################### # Removes unphysical low and high-frequency corrections def attenuation2(f, A, B, C): if min_freq < f < max_freq: return attenuation(f, A, B, C) else: return 0.0 # If freq > freq_limit no attenuation will be assumed ############################################################################## ##### Import Vt data ###### dataPath = '/home/pwb504/Documents/PhD/THz_TAW_2017/EMP_Analysis/Raw_Data_THz_2017/' shot1 = np.arange(3,15) shot2 = np.arange(16,59) shot3 = np.arange(60, 157) #shot4 = np.array([158]) shot5 = np.arange(161, 178) shots = np.concatenate([shot1, shot2, shot3, shot5]) shots=[5] channels = [2, 3, 4] # channel 3 and 4 are B-dots, so need to do different technique for C2 (Ddot) for c in range(len(channels)): for s in range(len(shots)): file = open(dataPath + 'ch' + str(channels[c]) + '_' + str(shots[s]) + '.csv') reader = csv.reader(file) # reads in file t = [] V = [] for row in reader: t.append(row[3]) # takes col 3 and col 4 from file and puts them into lists t and V V.append(row[4]) t = [float(j) for j in t] # change data to float format (prev. text) V = [float(k) for k in V] ####### CROP ############### crop_t = [] crop_V = [] # crop t between 380 and 800 ns for l in range(len(t)): indices = [] if 380E-9 < float(t[l]) < 800E-9: crop_t.append(float(t[l])) indices.append(l) for k in range(len(indices)): crop_V.append(float(V[indices[k]])) t = crop_t V = crop_V ######## Calculate Fourier Transform ############ signal = np.array(V, dtype=float) # convert signal to an array f_signal = np.fft.fft(signal) # compute frequencies for x-axis timestep = t[2] - t[1] n = signal.size freq = np.fft.fftfreq(n, d=timestep) ############################################################################### ############### Calculate attenuation and apply ############################### ################## Correct the FFT data ############################# attenuation_array = [] # use to check you get the right dB attenuation corrected = [] # Array of corrected data (i.e. unattenuated) for m in range(len(freq)): att_dB_100ft_new = attenuation2(abs(freq[m]), *popt) # attenuation per 100ft for freq[m] att_dB = (att_dB_100ft_new / 100) * (cable_length / 0.3048) # 0.3048 m in 1 foot (attenuation for freq[m]) V_in = f_signal[m] / (10**(att_dB / 20.0)) corrected.append(V_in) attenuation_array.append(att_dB) ### Remove low frequencies entirely - HIGH-PASS #### high_pass=[] for j in range(len(corrected)): if freq[j]min_freq: high_pass.append(corrected[j]) corrected = high_pass ##### Inverse Transform ###### V_corrected = np.fft.ifft(corrected) V_orig = np.fft.ifft(f_signal) V_corrected = V_corrected.real ############################################ ########## Find E-field #################### if channels[c] == 2: E_values = [0.0] # Assuming E0=0, have seed = E(t1) for i in range(1, len(t)): E_prev = E_values[i-1] EA = 0.5 * (V_corrected[i] + V_corrected[i-1]) * (t[i] - t[i-1]) + E_prev E_values.append(EA) E_values = [j / (50 * Ddot_Area * 8.85E-12) for j in E_values] ##### Energy in E-field ##### # Method 1: Integrate under square of E-field profile, multiply by const and multiply by coil area Flux = 0.0 for j in range(1, len(t)): E2_area = (0.5 * const / (9E16)) * ( (np.abs(E_values[j])**2) + (np.abs(E_values[j-1])**2) ) * (t[j] - t[j-1]) # Trapezium rule * const / c^2 Flux += E2_area # PF = instantaneous Poynting Flux (i.e PF at time t) EMP_Energy = Flux * Ddot_Area Energy_Array.append(EMP_Energy) ############################################ ########## Find B-field #################### else: ##### Correct for Balun 100G matching box ##### V_corrected2 = [m * (10**(8.0/20.0)) for m in V_corrected] # 8dB insertion loss through Balun ################# B_values = [0.0] # Assuming seed field is B0=0 for i in range(1, len(t)): B_prev = B_values[i-1] BA = 0.5 * (V_corrected2[i] + V_corrected2[i-1]) * (t[i] - t[i-1]) + B_prev # Add Balun correction B_values.append(BA) B_values = [j / Coil_Area for j in B_values] ########################################### ############## Energy ##################### ##### Energy in B-field ##### # Method 1: Integrate under square of B-field profile, multiply by const and multiply by coil area Flux1 = 0.0 for j in range(1, len(t)): B2_area = 0.5 * const * ( (np.abs(B_values[j])**2) + (np.abs(B_values[j-1])**2) ) * (t[j] - t[j-1]) # Trapezium rule * const Flux1 += B2_area # PF = instantaneous Poynting Flux (i.e PF at time t) EMP_Energy = Flux1 * Coil_Area Energy_Array.append(EMP_Energy) ### save to file for later ### # datafile path datafile_Path = '/home/pwb504/Documents/PhD/THz_TAW_2017/EMP_Analysis/Cable Corrected FFT Data/' + 'mod_ch' + str(channels[c]) + '_'+ str(shots[s]) +'.txt' # datafile id datafile_ID = open(datafile_Path, 'w+') # datafile is now open # Prepare data to write data_array2 = corrected # CORRECTED FFT data data_array1 = freq # frequency array # write np.savetxt(datafile_ID, np.transpose([data_array1, data_array2])) # transpose so saved as a column # can also achieve via freq.T datafile_ID.close() # close file ########### Save EMP Energy Data ################################################## # datafile path datafile_Path2 = '/home/pwb504/Documents/PhD/THz_TAW_2017/EMP_Analysis/EMP Energy/EMP_Energies.txt' # datafile id datafile_ID2 = open(datafile_Path2, 'w+') # datafile is now open # Prepare data to write Energy_Array_Ch2 = Energy_Array[0:len(shots)] # EMP Energies Energy_Array_Ch3 = Energy_Array[len(shots):2*len(shots)] Energy_Array_Ch4 = Energy_Array[2*len(shots):3*len(shots)] # write print len(shots), len(Energy_Array_Ch2), len(Energy_Array_Ch3) np.savetxt(datafile_ID2, np.transpose([shots, Energy_Array_Ch2, Energy_Array_Ch3, Energy_Array_Ch4])) # transpose so saved as a column datafile_ID2.close() # can also achieve via freq.T datafile_ID2.close() # close file