你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-08 07:24:19 +08:00
Updated the (main-like) function
这个提交包含在：
AdittyaPal
@@ -1,7 +1,9 @@
|
||||
|
||||
import numpy as np
|
||||
cimport cython
|
||||
from libc.math cimport floor, ceil, round, sqrt, tan, cos, sin, atan2, abs, pow, exp, M_PI
|
||||
from libc.stdio cimport FILE, fopen, fwrite, fclose
|
||||
from libc.string cimport strcmp
|
||||
from cython.parallel import prange
|
||||
|
||||
@cython.cdivision(True)
|
||||
@@ -101,6 +103,21 @@ cdef inline double getTanValue(int[:] integers, double[:] dr):
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef int[:, :] get_mZ(int m_x, int m_y, double theta, double[:] Delta_r):
|
||||
'''
|
||||
Method to get the arrays to perform a binary search to determine a rational number, m_z, closest to the real number, m_z, to get the desired tan Theta value.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
m_x and m_y, int : integers approximating the rational angle to the tan value of phi
|
||||
theta, float : polar angle of the incident plane wave (in radians) to be approximated to a rational angle
|
||||
Delta_r, float array : projections of the propagation vector along the different coordinate axes
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
_integers, int array : a two dimensional sequence of the three integers [m_x, m_y, m_z] to perform a binary search to determine the value of m_z within the given limits.
|
||||
'''
|
||||
cdef double m_z = 0
|
||||
m_z = sqrt((m_x/Delta_r[0])*(m_x/Delta_r[0]) + (m_y/Delta_r[1])*(m_y/Delta_r[1]))/(tan(theta)/Delta_r[2]) #get an estimate of the m_z value
|
||||
return np.array([[m_x, m_y, floor(m_z)],
|
||||
@@ -119,10 +136,10 @@ cdef int[:] getTheta(int m_x, int m_y, double theta, double Delta_theta, double[
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
m_x and m_y, int : integers approximating the rational angle to the tan value of phi
|
||||
theta, float : polar angle of the incident plane wave (in radians) to be approximated to a rational angle
|
||||
Delta_theta, float : permissible error in the rational angle approximation to theta (in radians)
|
||||
Delta_r, float : projections of the propagation vector along the different coordinate axes
|
||||
m_x and m_y, int : integers approximating the rational angle to the tan value of phi
|
||||
theta, float : polar angle of the incident plane wave (in radians) to be approximated to a rational angle
|
||||
Delta_theta, float : permissible error in the rational angle approximation to theta (in radians)
|
||||
Delta_r, float array : projections of the propagation vector along the different coordinate axes
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
@@ -169,7 +186,8 @@ def getIntegerForAngles(double phi, double Delta_phi, double theta, double Delta
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
integers, int array : sequence of the three integers [m_x, m_y, m_z]
|
||||
directions, int array : the integers specifying the direction of propagation of the plane wave along the three coordinate axes
|
||||
integers, int array : sequence of the three integers [m_x, m_y, m_z]
|
||||
'''
|
||||
cdef double required_ratio_phi, tolerance_phi = 0.0
|
||||
cdef int m_x, m_y, m_z = 0
|
||||
@@ -210,6 +228,25 @@ def getIntegerForAngles(double phi, double Delta_phi, double theta, double Delta
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, :, :, :, ::1] updateHFields(double[:] coefficients, int n_x, int n_y, int n_z,
|
||||
double[:, :, :, :, ::1] fields, int waveID):
|
||||
'''
|
||||
Method to update the magnetic fields in the TFSF grid uding the electric fields from the previous update.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
coefficients, float : stores the coefficients of the fields in the update equation for the magnetic field
|
||||
n_x, int : stores the number of grid cells along the x axis of the TFSF box
|
||||
n_y, int : stores the number of grid cells along the y axis of the TFSF box
|
||||
n_z, int : stores the number of grid cells along the z axis of the TFSF box
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
waveID, int : stores the index number of the field in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
fields, double array : returns the updated magnetic fields for the grid cells over the TFSF box at particular indices
|
||||
'''
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
|
||||
cdef double c0 = coefficients[0]
|
||||
@@ -245,6 +282,25 @@ cdef double[:, :, :, :, ::1] updateHFields(double[:] coefficients, int n_x, int
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, :, :, :, ::1] updateEFields(double[:] coefficients, int n_x, int n_y, int n_z,
|
||||
double[:, :, :, :, ::1] fields, int waveID):
|
||||
'''
|
||||
Method to update the electric fields in the TFSF grid uding the electric fields from the previous update.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
coefficients, float : stores the coefficients of the fields in the update equation for the electric field
|
||||
n_x, int : stores the number of grid cells along the x axis of the TFSF box
|
||||
n_y, int : stores the number of grid cells along the y axis of the TFSF box
|
||||
n_z, int : stores the number of grid cells along the z axis of the TFSF box
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
waveID, int : stores the index number of the field in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
fields, double array : returns the updated electric fields for the grid cells over the TFSF box at particular indices
|
||||
'''
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
|
||||
cdef double c0 = coefficients[0]
|
||||
@@ -285,6 +341,27 @@ cdef double[:, :, :, :, ::1] updateEFields(double[:] coefficients, int n_x, int
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, :, :, ::1] implementABC(double[:, :, :, ::1] face_fields, double abccoef,
|
||||
int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int waveID):
|
||||
'''
|
||||
Method to implement the absorbing boundary conditions to the FDTD grid.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
face_fields, double array : stores the fields on the six faces of the FDTD Box
|
||||
abccoef, double : stores the absorbing boundary coefficient so that there are no reflection from the edge of the grid
|
||||
n_x, int : stores the number of grid cells along the x axis of the TFSF box
|
||||
n_y, int : stores the number of grid cells along the y axis of the TFSF box
|
||||
n_z, int : stores the number of grid cells along the z axis of the TFSF box
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
waveID, int : stores the index number of the field in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
face_fields, double array : returns the updated fields on the surfaces of the FDTD Box
|
||||
|
||||
'''
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
|
||||
cdef double[:, :, ::1] e_x = fields[waveID, 0, :, :, :]
|
||||
@@ -357,6 +434,27 @@ cdef double[:, :, :, ::1] implementABC(double[:, :, :, ::1] face_fields, double
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, :, :, :, ::1] applyTFSFMagnetic(double[:] coefficients, double[:, :] e1D, int[:] m,
|
||||
double[:, :, :, :, ::1]fields, int[:, :] corners, int waveID):
|
||||
'''
|
||||
Method to implement the total field-scattered field formulation for the magnetic field on the edge of the TF/SF region of the TFSF Box.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
coefficients, float : stores the coefficients of the fields in the TFSF assignment equation for the magnetic field
|
||||
e1D, double array : array to store the electric fields of the one dimensional representation of the plane wave in a direction along which the wave propagates
|
||||
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles, for assignment of the correct element
|
||||
to the three dimensional FDTD grid from the one dimensional representation, last element stores max(m_x, m_y, m_z)
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries
|
||||
waveID, int : stores the index number of the field in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
fields, double array : returns the updated fields atre applying the TF/Sf formulation at the boundary of the TF/SF surface
|
||||
|
||||
'''
|
||||
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
cdef int x_start = corners[0, 0]
|
||||
@@ -462,6 +560,27 @@ cdef double[:, :, :, :, ::1] applyTFSFMagnetic(double[:] coefficients, double[:,
|
||||
cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:, :] h1D, int[:] m,
|
||||
double[:, :, :, :, ::1] fields, int[:, :] corners, int waveID):
|
||||
|
||||
'''
|
||||
Method to implement the total field-scattered field formulation for the electric field on the edge of the TF/SF region of the TFSF Box.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
coefficients, float : stores the coefficients of the fields in the TFSF assignment equation for the electric field
|
||||
h1D, double array : array to store the magnetic fields of the one dimensional representation of the plane wave in a direction along which the wave propagates
|
||||
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles, for assignment of the correct element
|
||||
to the three dimensional FDTD grid from the one dimensional representation, last element stores max(m_x, m_y, m_z)
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries
|
||||
waveID, int : stores the index number of the field in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
fields, double array : returns the updated fields atre applying the TF/Sf formulation at the boundary of the TF/SF surface
|
||||
|
||||
'''
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
cdef int x_start = corners[0, 0]
|
||||
cdef int y_start = corners[0, 1]
|
||||
@@ -561,7 +680,13 @@ cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:,
|
||||
e_x[i, j, k] -= coefficients[2] * h1D[1, index]
|
||||
return fields
|
||||
|
||||
'''
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, ::1] updateMagneticFields(int n, double[:] H_coefficients, double[:, ::1] H_fields,
|
||||
double[:, ::1] E_fields, int dimensions, int[:] m):
|
||||
'''
|
||||
Method to update the magnetic fields for the next time step using
|
||||
Equation 8 of DOI: 10.1109/LAWP.2009.2016851
|
||||
__________________________
|
||||
@@ -572,20 +697,17 @@ cdef double[:, :, :, :, ::1] applyTFSFElectric(double[:] coefficients, double[:,
|
||||
H_coefficients, double array : stores the coefficients of the fields in the update equation for the magnetic field
|
||||
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
|
||||
E_fields, double array : stores the electric fields of the DPW till temporal index time
|
||||
time, int : time index storing the current axis number which would be updated for the H_fields
|
||||
dimensions, int : stores the number of dimensions in which the simulation is run
|
||||
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles, for assignment of the correct element
|
||||
to the three dimensional FDTD grid from the one dimensional representation, last element stores max(m_x, m_y, m_z)
|
||||
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
H_fields, double array : magnetic field array with the axis entry for the current time added
|
||||
|
||||
'''
|
||||
@cython.cdivision(True)
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, ::1] updateMagneticFields(int n, double[:] H_coefficients, double[:, ::1] H_fields,
|
||||
double[:, ::1] E_fields, int dimensions, int[:] m):
|
||||
|
||||
'''
|
||||
cdef Py_ssize_t i, j = 0
|
||||
cdef int dim_mod1, dim_mod2 = 0
|
||||
for i in prange(dimensions, nogil=True, schedule='static'): #loop to update each component of the magnetic field
|
||||
@@ -613,7 +735,10 @@ cdef double[:, ::1] updateElectricFields(int n, double[:] E_coefficients, double
|
||||
E_coefficients, double array : stores the coefficients of the fields in the update equation for the electric field
|
||||
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
|
||||
E_fields, double array : stores the electric fields of the DPW till temporal index time
|
||||
time, int : time index storing the current axis number which would be updated for the H_fields
|
||||
dimensions, int : stores the number of dimensions in which the simulation is run
|
||||
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles, for assignment of the correct element
|
||||
to the three dimensional FDTD grid from the one dimensional representation, last element stores max(m_x, m_y, m_z)
|
||||
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
@@ -633,29 +758,81 @@ cdef double[:, ::1] updateElectricFields(int n, double[:] E_coefficients, double
|
||||
return E_fields
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef inline double getSource(double time, double ppw):
|
||||
cdef double getSource(double time, double ppw, char* wavetype, double dt):
|
||||
'''
|
||||
Method to get the magnitude of the source field in the direction perpendicular to the propagation of the plane wave
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
time, float : time at which the magnitude of the source is calculated
|
||||
ppw, int : points per wavelength for the wave source
|
||||
time, float : time at which the magnitude of the source is calculated
|
||||
ppw, int : points per wavelength for the wave source
|
||||
wavetype, string : stores the type of waveform whose magnitude should be returned
|
||||
dt, double : the time upto which the wave should exist in a impulse delta pulse
|
||||
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
sourceMagnitude, double array : magnitude of the source for the requested indices at the current time
|
||||
sourceMagnitude, double : magnitude of the source for the requested indices at the current time
|
||||
'''
|
||||
return (1.0-2.0*(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0))
|
||||
)*exp(-(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0))) # define a Ricker wave source
|
||||
# Waveforms
|
||||
if (strcmp(wavetype, "gaussian") == 0):
|
||||
return exp(-(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandot") == 0 or strcmp(wavetype, "gaussianprime") == 0):
|
||||
return 4 * M_PI/ppw * (M_PI*(time/<double>ppw - 1.0)
|
||||
) * exp(-(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotnorm") == 0):
|
||||
return 2 * (M_PI*(time/<double>ppw - 1.0)
|
||||
) * exp(-(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotdot") == 0 or strcmp(wavetype, "gaussiandoubleprime") == 0):
|
||||
return 2 * M_PI*M_PI/(ppw*ppw) * (1.0-2.0*(M_PI*(time/ppw - 1.0))*(M_PI*(time/ppw - 1.0))
|
||||
) * exp(-(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotdotnorm") == 0 or strcmp(wavetype, "ricker") == 0):
|
||||
return (1.0 - 2.0*(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0))
|
||||
) * exp(-(M_PI*(time/<double>ppw - 1.0)) * (M_PI*(time/<double>ppw - 1.0))) # define a Ricker wave source
|
||||
|
||||
elif (strcmp(wavetype, "sine") == 0):
|
||||
return sin(2 * M_PI * time/<double>ppw)
|
||||
|
||||
elif (strcmp(wavetype, "contsine") == 0):
|
||||
return min(0.25 * time/<double>ppw, 1) * sin(2 * M_PI * time/<double>ppw)
|
||||
|
||||
elif (strcmp(wavetype, "impulse") == 0):
|
||||
if (time < dt): # time < dt condition required to do impulsive magnetic dipole
|
||||
return 1
|
||||
else:
|
||||
return 0
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef double[:, :, :, :, ::1] superImposeFields(planeWaves, int n_x, int n_y, int n_z, int dimensions,
|
||||
double[:, :, :, :, ::1] fields, int iterate):
|
||||
'''
|
||||
Method to superimpose multiple plane waves if the user requires multiple plane wave incidences.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
planeWaves, list : list of class DiscretePlaneWave() objects, each representing a particular plane wave
|
||||
n_x, int : stores the number of grid cells along the x axis of the TFSF box
|
||||
n_y, int : stores the number of grid cells along the y axis of the TFSF box
|
||||
n_z, int : stores the number of grid cells along the z axis of the TFSF box
|
||||
dimensions, int : stores the number of dimensions in which the simulation is run
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
iterate, int : stores the number of incident plane waves in case of multiple plane waves
|
||||
__________________________
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
fields, double array : returns the superimposed fields in the TFSF box as the fields in the last index.
|
||||
|
||||
'''
|
||||
cdef Py_ssize_t i, j, k, l, n = 0
|
||||
cdef int x_start, y_start, z_start
|
||||
for n in range(iterate):
|
||||
@@ -675,45 +852,71 @@ cdef double[:, :, :, :, ::1] superImposeFields(planeWaves, int n_x, int n_y, int
|
||||
fields[iterate, l, i, j, k] += fields[n, l, abs(x_start-i), abs(y_start-j), abs(z_start-k)]
|
||||
|
||||
return fields
|
||||
'''
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef void takeSnapshot3D(double[:, :, ::1] field, char* filename):
|
||||
cdef FILE *fptr = fopen(filename, "wb")
|
||||
cdef Py_ssize_t n_x = field.shape[0]
|
||||
cdef double[:, ::1] slice_view
|
||||
cdef Py_ssize_t i = 0
|
||||
if fptr is NULL:
|
||||
raise ValueError("Failed to open the file.")
|
||||
else:
|
||||
for i in range(n_x):
|
||||
slice_view = field[i, ...]
|
||||
fwrite(&slice_view[0, 0], sizeof(double), slice_view.size, fptr)
|
||||
fclose(fptr)
|
||||
'''
|
||||
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
@cython.cdivision(True)
|
||||
cdef void takeSnapshot3D(double[:, :, ::1] field, char* filename):
|
||||
'''
|
||||
Method to write the fields of the plane wave simulation at a particular time step.
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices of the TFSF box at a particular time step
|
||||
filename, string : specifies the file location where the fields are to be written
|
||||
|
||||
'''
|
||||
cdef FILE *fptr = fopen(filename, "wb")
|
||||
fwrite(&field[0, 0, 0], sizeof(double), field.size, fptr)
|
||||
fclose(fptr)
|
||||
|
||||
|
||||
def getGridFields(planeWave, double[:] C, double[:] D, int snapshot,
|
||||
int n_x, int n_y, int n_z, double[:, :, :, :, ::1] fields, int[:, ::1] corners,
|
||||
int time_duration, double[:, :, :, ::1] face_fields,
|
||||
double abccoef, double dt, int noOfWaves, double c, double ppw, int dimensions):
|
||||
|
||||
double abccoef, double dt, int noOfWaves, double c, double ppw,
|
||||
int dimensions, filepath, snapFieldNumber, snapFieldSteps):
|
||||
'''
|
||||
Method to implement the time stepping loop of the FDTD simulation
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
planeWave, list : list of class DiscretePlaneWave() objects, each representing a particular plane wave
|
||||
C, double array : stores the coefficients of the fields for the update equation of the electric fields
|
||||
D, double array : stores the coefficients of the fields for the update equation of the magnetic fields
|
||||
snapshot, int : stores the interval at which a snapshot of the fields is stored
|
||||
n_x, int : stores the number of grid cells along the x axis of the TFSF box
|
||||
n_y, int : stores the number of grid cells along the y axis of the TFSF box
|
||||
n_z, int : stores the number of grid cells along the z axis of the TFSF box
|
||||
fields, double array : stores the fields for the grid cells over the TFSF box at particular indices in the order
|
||||
E_x, E_y, E_z, H_x, H_y, H_z
|
||||
corners, int array : stores the coordinates of the cornets of the total field/scattered field boundaries
|
||||
time_dimension, int : stores the number of time steps over which the FDTD simulation is run
|
||||
face_fields, double array : stores the fields on the six faces of the FDTD Box
|
||||
abccoef, double : stores the absorbing boundary coefficient so that there are no reflection from the edge of the grid
|
||||
noOfWaves, int : stores the number of plane waves to be incident on the TFSF Box
|
||||
c, double : speed of light
|
||||
ppw, double : stores the number of points per wavelength for the requested source
|
||||
dimensions, int : stores the number of dimensions in which the simulation is run
|
||||
filepath, string : specifies the file location where the fields are to be written
|
||||
snapFieldNumber, list: list to store which fields are to be written in the file
|
||||
snapFieldSteps, list : list to store any specific time steps at which snapshots are to be taken.
|
||||
|
||||
'''
|
||||
cdef double real_time = 0.0 # initialize the real time to zero
|
||||
# Py_ssize_t is the proper C type for Python array indices
|
||||
cdef Py_ssize_t i, j, r = 0 #E and H are the one-dimensional waves
|
||||
cdef Py_ssize_t i, j, k, l, m = 0 #E and H are the one-dimensional waves
|
||||
for i in range(time_duration):
|
||||
for j in range(noOfWaves):
|
||||
for k in range(dimensions):
|
||||
for l in range(planeWave[j].m[dimensions]): #loop to assign the source values of magnetic field to the first few gridpoints
|
||||
planeWave[j].H_fields[k, l] = planeWave[j].projections[k]*getSource(real_time - (l+(
|
||||
planeWave[j].m[(k+1)%dimensions]+planeWave[j].m[(k+2)%dimensions]
|
||||
)*0.5)*planeWave[j].ds/c, ppw)
|
||||
)*0.5)*planeWave[j].ds/c, ppw,
|
||||
planeWave[j].waveformID, planeWave[j].dt)
|
||||
|
||||
planeWave[j].H_fields = updateMagneticFields(planeWave[j].length, D, planeWave[j].H_fields,
|
||||
planeWave[j].E_fields, dimensions, planeWave[j].m) #update the magnetic fields in the rest of the DPW grid
|
||||
@@ -730,8 +933,10 @@ def getGridFields(planeWave, double[:] C, double[:] D, int snapshot,
|
||||
|
||||
real_time += dt #take a half time step because the magnetic field and electric field grids are staggered in time
|
||||
|
||||
if (i % snapshot == 0):
|
||||
if (i % snapshot == 0 or i in snapFieldSteps):
|
||||
fields = superImposeFields(planeWave, n_x, n_y, n_z, dimensions, np.asarray(fields), noOfWaves)
|
||||
takeSnapshot3D(fields[noOfWaves, 0, :, :, :], f'./snapshots/electric_2{i}.dat'.encode('UTF-8'))
|
||||
fields[noOfWaves, :, :, :, :] = 0
|
||||
for m in snapFieldNumber:
|
||||
takeSnapshot3D(fields[noOfWaves, m, :, :, :], (f'{filepath}_{m}{i}.dat').encode('UTF-8'))
|
||||
|
||||
fields[noOfWaves, :, :, :, :] = 0 #erase the cells bearing the last index so that it does not interfere with the superimposiion fields at some other time instance
|
||||
|
||||
|
||||
在新工单中引用
屏蔽一个用户