修改了两处不严谨的地方，gfortran环境下会warning，1：get_non_uniformgrid.f90中，IF判断中IF(SourceLength/GridSize.gt.51)改SourceLength/GridSize为整型SourceGridNum；2：替换旧写法PAUSE语句为READ(*,*)

Merge branch 'master' of https://git.em3d.cn/openSource/tem3dfdtd-open
修正一处bug，VTK_fortran.f90 格式化输出格式不规范，无法在gfortran环境中编译，现已改为标准格式
2026-08-18 18:48:15 +08:00 · 2026-08-18 16:21:44 +08:00 · 2026-08-18 16:19:47 +08:00 · 2026-08-17 21:00:41 +08:00 · 2026-08-14 20:54:09 +08:00 · 2026-08-14 20:51:58 +08:00
@@ -252,3 +252,4 @@ ModelManifest.xml
 /tem3dfdtd/HzCoordinate.dat
 /tem3dfdtd/Ground-Line=081.dat
 tem3dfdtd/fort.5141
+.vscode/settings.json
@@ -1,4 +1,545 @@
-# tem3dfdtd-open
+# tem3dfdtd

-#### Description
-Please visit https://em3d.cn for more information.
+[中文](./README.md) | [English](./README.en.md)
+
+[!NOTE]
+
+> All contents of this document are automatically translated by AI without human verification. Please refer to the original Chinese version if you need accurate information.
+
+# 3D Transient Electromagnetic FDTD Forward Modeling Program (tem3dfdtd)
+
+This program performs forward modeling of three-dimensional transient electromagnetic (TEM) responses based on the **FDTD (finite-difference time-domain) method**. The core algorithm adopts the Du Fort–Frankel method improved by Wang–Hohmann (1993): it iteratively solves the magnetic field H on a Yee grid, and ensures the time stability of the explicit iteration by introducing a **fictitious permittivity**. It uses **conformal grid techniques** to handle undulating terrain and arbitrarily shaped anomalous bodies (described by surface triangular meshes). It supports a rectangular loop source and can compute ground TEM, semi-airborne (SATEM), and airborne (ATEM) modes.
+
+The overall framework and the three core techniques are derived from the following works:
+
+- **Overall framework and core iterative algorithm** (loop-source TEM three-dimensional finite-difference time-domain forward modeling considering the turn-off time, the Wang–Hohmann improved Du Fort–Frankel method, the fictitious permittivity, and the source waveform including the turn-off time); for the theoretical details, see reference [1];
+- **CPML absorbing boundary** (the CPML absorbing boundary for the low-frequency-approximated Maxwell equations in transient electromagnetics and its application method); for the theoretical details, see reference [2];
+- **Conformal grid technique** (introducing arbitrarily complex-shaped structures into the Yee grid computation via a ray-tracing method); for the theoretical details, see reference [3].
+
+The specific principles, formula derivations, and implementation details of each method are given in references [1]–[3]. The authors of these papers are all people who have made outstanding contributions to this open-source project.
+
+Code structure: `main.f90` (main program), `module/` (global parameters and modules), `lib/` (functional subroutines).
+
+---
+
+## Table of Contents
+
+- [tem3dfdtd](#tem3dfdtd)
+- [3D Transient Electromagnetic FDTD Forward Modeling Program (tem3dfdtd)](#3d-transient-electromagnetic-fdtd-forward-modeling-program-tem3dfdtd)
+  - [Table of Contents](#table-of-contents)
+  - [1. Environment and Compilation](#1-environment-and-compilation)
+    - [1.1 Environment Requirements](#11-environment-requirements)
+    - [1.2 Build and Run with VS2019 (Recommended)](#12-build-and-run-with-vs2019-recommended)
+    - [1.3 Command-Line Compilation (Optional, Not Recommended on Windows)](#13-command-line-compilation-optional-not-recommended-on-windows)
+    - [1.4 Linux Environment Compilation (makefile)](#14-linux-environment-compilation-makefile)
+  - [2. Program Flow and Module Structure](#2-program-flow-and-module-structure)
+    - [Module 1: Program Control and Parameter Input](#module-1-program-control-and-parameter-input)
+    - [Module 2: Yee Grid Generation](#module-2-yee-grid-generation)
+    - [Module 3: Electrical Parameter Construction](#module-3-electrical-parameter-construction)
+    - [Module 4: Excitation Source and Time Series](#module-4-excitation-source-and-time-series)
+    - [Module 5: FDTD Electromagnetic Field Computation](#module-5-fdtd-electromagnetic-field-computation)
+    - [Module 6: Output](#module-6-output)
+  - [3. Input File Formats](#3-input-file-formats)
+    - [3.1 Parameter Control File input.dat](#31-parameter-control-file-inputdat)
+    - [3.2 Terrain Mesh File Complex_Terrain](#32-terrain-mesh-file-complex_terrain)
+    - [3.3 Anomalous Body Mesh File Complex_anomalous](#33-anomalous-body-mesh-file-complex_anomalous)
+    - [3.4 Modeling Notes](#34-modeling-notes)
+  - [4. CPML Absorbing Boundary](#4-cpml-absorbing-boundary)
+    - [4.1 Implementation Location](#41-implementation-location)
+    - [4.2 Parameters and Their Meanings](#42-parameters-and-their-meanings)
+    - [4.3 Usage Notes](#43-usage-notes)
+    - [4.4 Comparison with the Dirichlet Boundary](#44-comparison-with-the-dirichlet-boundary)
+  - [5. Output File Description](#5-output-file-description)
+  - [6. Quick Decay-Curve Plotting (TEM_decay_plot.py)](#6-quick-decay-curve-plotting-tem_decay_plotpy)
+    - [6.1 Usage](#61-usage)
+    - [6.2 Main Adjustable Parameters (the "User parameters" section at the top of the script)](#62-main-adjustable-parameters-the-user-parameters-section-at-the-top-of-the-script)
+  - [7. Quick Start (the example bundled in this directory)](#7-quick-start-the-example-bundled-in-this-directory)
+  - [8. FAQ](#8-faq)
+  - [9. References](#9-references)
+  - [10. Contributors](#10-contributors)
+  - [11. Statement](#11-statement)
+
+---
+
+## 1. Environment and Compilation
+
+### 1.1 Environment Requirements
+
+| Item             | Requirement                                                                                                                      |
+| ---------------- | -------------------------------------------------------------------------------------------------------------------------------- |
+| Operating system | Windows 10/11 (64-bit) or Linux (x86_64); the domestic operating system deepin 25 is especially recommended                      |
+| IDE              | Windows: Visual Studio 2019 or later; Linux/deepin: make;                                                                        |
+| Fortran compiler | Intel oneAPI Fortran (ifort/ifx);                                                                                                |
+| Parallel support | OpenMP (multi-core CPU acceleration); additional GPU support is provided in the commercial version, please visit https://em3d.cn |
+
+### 1.2 Build and Run with VS2019 (Recommended)
+
+1. Install Visual Studio 2019 (select the "Desktop development with C++" workload) and Intel oneAPI (select the "Visual Studio integration for Intel Fortran Compiler" during installation).
+2. Double-click to open the project file **`tem3dfdtd.sln`** (an Intel Fortran project, corresponding to
+   `tem3dfdtd.vfproj`). The project already includes all source files. Configuration notes:
+   - `Debug | x64` / `Release | x64`: use the **ifx** compiler (oneAPI default), **recommended**;
+   - `Debug | Win32` / `Release | Win32`: use the ifort compiler.
+3. Select the **`Release | x64`** configuration and click **Build → Build Solution**.
+4. Preparation before running: the program looks for `input.dat` and the mesh files in the working directory, so place
+   `input.dat`, `Complex_Terrain.dat/.stl`, and `Complex_anomalous.dat/.stl` in the
+   **tem3dfdtd-open\tem3dfdtd** directory (or specify it via "Project Properties → Debugging → Working Directory").
+5. Run `tem3dfdtd.exe` directly, or press F5 to debug-run.
+
+> **Note**: the program uses OpenMP and requires Intel's runtime library `libiomp5md.dll` at run time
+> (located in the Intel oneAPI installation directory `bin/`). If you get a missing-DLL error, copy it next to
+> the exe (a copy is already placed in this directory), or add its directory to the system PATH.
+
+### 1.3 Command-Line Compilation (Optional, Not Recommended on Windows)
+
+In an "Intel oneAPI Command Prompt" environment:
+
+```bat
+ifx -c -O2 -Qopenmp module\*.f90
+ifx -c -O2 -Qopenmp -Qopenmp lib\*.f90 main.f90
+ifx -O2 -Qopenmp *.obj -o tem3dfdtd.exe
+```
+
+(Linking requires MSVC's link.exe and the Windows SDK libraries; it is recommended to use VS's "Developer Command Prompt"
+and add Intel oneAPI's `bin` directory to PATH.)
+
+### 1.4 Linux Environment Compilation (makefile)
+
+To compile under Linux you need to write a `makefile`. Please write the makefile yourself according to your operating system's configuration and requirements, and compile with `make` (using the makefile).
+
+**Notes and cautions**
+
+- Differences from the Windows version: under Linux the executable name is custom; if not specified it defaults to `main.exe`. The working directory must likewise contain
+  `input.dat` and the `Complex_Terrain.*`, `Complex_anomalous.*` mesh files.
+
+---
+
+## 2. Program Flow and Module Structure
+
+Main program execution order (`main.f90`):
+
+```
+GETDATA → CHECKPARAMETERS → MEMORY_USE_ESTIMATION → ALLOCATEMEMORY
+→ GET_NON_UNIFORMGRID → ZERO → GET_COORDINATES
+→ Get_Receiver_Gridlabel → RES_CONFIGURE → TIME_SERIOUS
+→ Get_eps_r → (if Logic_PML=1) Get_pml_parameters → Get_mstop
+→ GetSourcePosition → Iteration → FREE_MEMORY
+```
+
+By function, it is divided into the following 6 modules:
+
+### Module 1: Program Control and Parameter Input
+
+| File                            | Function                                                                                                                             |
+| ------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------ |
+| `main.f90`                      | Main program, controls the whole computation flow                                                                                    |
+| `lib/getdata.f90`               | Reads the parameter control file `input.dat`; detects the existence of terrain/anomalous-body mesh files and selects the read format |
+| `lib/checkparameters.f90`       | Echoes the read-in computation parameters to `logfile.log` for manual checking                                                       |
+| `lib/memory-use-estimation.f90` | Estimates the required memory based on the grid size and prints a hint                                                               |
+
+### Module 2: Yee Grid Generation
+
+| File                          | Function                                                                                                                                       |
+| ----------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------- |
+| `lib/allocatememory.f90`      | Dynamically allocates all global arrays according to the input parameters (including the CPML memory-variable arrays, only when `Logic_PML=1`) |
+| `lib/get_non_uniformgrid.f90` | Generates the non-uniform grid in the x/y/z directions (uniform in the core area + expanded outward by a factor of 1.3)                        |
+| `lib/get_coordinates.f90`     | Computes the coordinates of all grid nodes (including Yee nodes); coordinates take the source center as origin                                 |
+| `lib/zero.f90`                | Initializes all electromagnetic field arrays to 0; sets den\_\* to 1, c_h_zz to 0, and clears the CPML memory variables                        |
+
+### Module 3: Electrical Parameter Construction
+
+| File                                | Function                                                                                                                                                                                                                                                                                        |
+| ----------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `lib/resistivity-configuration.f90` | Builds the model conductivity: without terrain, assigns background conductivity + block-shaped anomalous bodies; with terrain, calls the conformal mesh; finally distributes the conductivity onto edges in the x/y/z directions and writes `conductivity.vtk`                                  |
+| `lib/Terrain_conformal.f90`         | **Terrain conformal mesh**: reads the terrain triangular mesh from `Complex_Terrain.dat/.stl`, fills the equivalent conductivity of each edge along the x/y/z directions using ray–triangle intersection (Möller–Trumbore algorithm), handling undulating terrain and the air/stratum interface |
+| `lib/Anomalous_conformal.f90`       | **Anomalous-body conformal mesh**: reads the anomalous-body surface mesh from `Complex_anomalous.dat/.stl`, and fills the anomalous-body conductivity (`tao_abnormal`) onto edges using the same ray-intersection method                                                                        |
+
+> **Note**: when a terrain file exists, in the terrain branch the air conductivity is `AIR_CONDUCTIVITY = 1e-6 S/m`,
+> and the stratum conductivity is the **2nd conductivity in the block-shaped anomalous body parameters** `TAR_CONDUCTIVITY(2)` (see the second group of data from line 9 of input.dat). In this case the block-shaped anomalous bodies in input.dat do not take effect directly; instead the anomalous body is described by the `Complex_anomalous`
+> mesh file and assigned its conductivity via `tao_abnormal`.
+
+### Module 4: Excitation Source and Time Series
+
+| File                          | Function                                                                                                                                                                                        |
+| ----------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `lib/time-serious.f90`        | Generates the time series for the whole computation (including the source waveform), and corrects the iteration step count `NSTOP` according to `MAX_OFF_TIME`; writes `CTIME_TIXING_UPCOS.DAT` |
+| `lib/tixing-source-upcos.f90` | Source waveform: **trapezoidal + cosine-ramp** turn-off current waveform (commonly used, `SOURCE_TYPE = 'TIXING_UPCOS'`)                                                                        |
+| `lib/tixing-source.f90`       | Pure trapezoidal waveform source (`TIXING_RAMP`)                                                                                                                                                |
+| `lib/sin-source.f90`          | Half-sine waveform source (`HALF_SIN`)                                                                                                                                                          |
+| `lib/triangle-source.f90`     | Triangular waveform source (`TRIANGLE`)                                                                                                                                                         |
+| `lib/get-eps-r.f90`           | Computes the fictitious permittivity `EPS_R = 3·(Δt/Δx)²/μ₀` and the iteration coefficients, ensuring the stability of the explicit FDTD                                                        |
+| `lib/get-mstop.f90`           | Splits the total iteration into several "computation segments", each with its own cache allocation, for memory management                                                                       |
+
+### Module 5: FDTD Electromagnetic Field Computation
+
+| File                        | Function                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
+| --------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `lib/GetSourcePosition.f90` | Determines the loop source's position in the grid and marks the source edges (the current assignment region)                                                                                                                                                                                                                                                                                                                                                                                                                                 |
+| `lib/Iteration.f90`         | **Core iteration subroutine**: advances the time steps segment by segment, updates the Ex/Ey/Ez and Hz fields, and loads current according to the source waveform; at the end of each segment it computes Hz for each receiver point (obtained by weighted differencing of the Ex, Ey of the 8 nodes surrounding that point) and writes it to the result file. The boundary condition switches via `Logic_PML`: `1` adds the CPML memory-variable correction in each field update, `0` restores the original Dirichlet (zero-field) boundary |
+| `module/pml-parameters.f90` | **CPML absorbing boundary module** (Roden–Gedney convolutional PML): declares the σ/α/κ polynomial scaling parameters (ma=3, mb=1), the 26 memory-variable ψ arrays, the b*e/c_e convolution-coefficient arrays, and the den*\*(=1/κ) scaling arrays                                                                                                                                                                                                                                                                                         |
+| `lib/get-pml-paramters.f90` | Builds the σ/α/κ distributions on the six x/y/z boundary faces (polynomial decay from the boundary inward, E/H staggered sampling) and the den\_\* scaling arrays; called by main only when `Logic_PML=1`                                                                                                                                                                                                                                                                                                                                    |
+
+### Module 6: Output
+
+| File                                              | Function                                                                          |
+| ------------------------------------------------- | --------------------------------------------------------------------------------- |
+| `lib/Iteration.f90` (output part)                 | Writes the response files `dBzdt_1.txt`, `dBzdt_2.txt`, … for each receiver point |
+| `lib/resistivity-configuration.f90` (output part) | Writes the model conductivity distribution `conductivity.vtk`                     |
+| `lib/time-serious.f90` (output part)              | Writes the time series `CTIME_TIXING_UPCOS.DAT`                                   |
+| `lib/free-memory.f90`                             | Releases all dynamic memory after computation                                     |
+
+---
+
+## 3. Input File Formats
+
+The program requires the following files (all placed in the exe's working directory):
+
+| File                              | Required?             | Description                                                                                     |
+| --------------------------------- | --------------------- | ----------------------------------------------------------------------------------------------- |
+| `input.dat`                       | **Required**          | Computation parameter control file                                                              |
+| `Complex_anomalous.dat` or `.stl` | Optional (either one) | Anomalous-body surface triangular mesh; if absent, the model is treated as a uniform background |
+| `Complex_Terrain.dat` or `.stl`   | Optional (either one) | Terrain surface triangular mesh; if absent, no undulating terrain is considered                 |
+
+### 3.1 Parameter Control File input.dat
+
+**Free-format reading**, read line by line in order; a value may be followed by a `!` comment (whole-line or end-of-line comments). The following explains each line using the `input.dat` bundled in this directory as an example:
+
+| Line  | Example        | Meaning                                                                                                                                   |
+| ----- | -------------- | ----------------------------------------------------------------------------------------------------------------------------------------- |
+| 1     | `1`            | Computation mode `CAL_TYPE`: `1` = ground TEM, `2` = semi-airborne (SATEM)                                                                |
+| 2     | `500`          | Transmitting loop side length `SourceLength` (m)                                                                                          |
+| 3     | `101,101,100`  | Number of grids in the x, y, z directions `NX,NY,NZ`                                                                                      |
+| 4     | `1`            | **Boundary condition switch `Logic_PML`**: `1` = CPML absorbing boundary, `0` = original non-uniform-grid Dirichlet (zero-field) boundary |
+| 5     | `10,10,10`     | **PML layer count `PML_X,PML_Y,PML_Z`** (x, y, z directions; effective only when the switch = 1; ≥ 5 layers recommended)                  |
+| 6     | `25,25`        | x direction: core uniform-grid start/end interval indices `UniGridNumX1,UniGridNumX2`                                                     |
+| 7     | `25,25`        | y direction: core uniform-grid interval indices `UniGridNumY1,UniGridNumY2`                                                               |
+| 8     | `20,30`        | z direction: core uniform-grid interval indices `UniGridNumZ1,UniGridNumZ2`                                                               |
+| 9     | `20`           | Core-area uniform grid size `GridSize` (m)                                                                                                |
+| 10    | `0.01`         | Background medium conductivity `BACKGROUND_CONDUCTIVITY` (S/m)                                                                            |
+| 11    | `2`            | Number of block-shaped anomalous bodies `TEMP_II` (filled as prisms when no terrain; set 0 for a uniform model)                           |
+| 12–15 | see below      | **Parameters of the 1st block-shaped anomalous body**, 4 lines in total                                                                   |
+| 16–19 | see below      | **Parameters of the 2nd block-shaped anomalous body**, 4 lines in total                                                                   |
+| 20    | `4000000`      | Maximum number of iterations `NSTOP`                                                                                                      |
+| 21    | `90.101`       | Maximum computation time `MAX_OFF_TIME` (**unit: ms**)                                                                                    |
+| 22    | `1e-6,1e-9`    | Rising-edge duration and time step `RAISETIME, RAISESTEP` (s)                                                                             |
+| 23    | `60000e-6`     | Plateau-stage duration `WAVE` (s, i.e. 60 ms)                                                                                             |
+| 24    | `1e-7,1e-9`    | Falling-edge duration and time step `RAMP, RAMPSTEP` (s)                                                                                  |
+| 25    | `1e-9`         | Initial time step `TIMESTEP` (s)                                                                                                          |
+| 26    | `1`            | Transmitting current amplitude `AMP` (A)                                                                                                  |
+| 27    | `4.0`          | Anomalous-body conductivity `tao_abnormal` (S/m, used with the `Complex_anomalous` file)                                                  |
+| 28    | `TIXING_UPCOS` | Source type `SOURCE_TYPE`: `TIXING_UPCOS` / `TIXING_RAMP` / `HALF_SIN` / `TRIANGLE`                                                       |
+| 29    | `1`            | Number of receiver points `Point_Num`                                                                                                     |
+| 30    | `1`            | Index of the 1st receiver point                                                                                                           |
+| 31    | `0,0,0`        | Coordinates of the 1st receiver point (relative to the source center, m)                                                                  |
+| 32–33 | `2` / `x,y,z`  | (if there are extra measurement points) the 2nd receiver point (index + coordinates)                                                      |
+
+Each **block-shaped anomalous body** consists of 4 consecutive lines:
+
+| Line | Example | Meaning                                             |
+| ---- | ------- | --------------------------------------------------- |
+| ①    | `1,101` | x-direction grid start/end indices `TAR_X1, TAR_X2` |
+| ②    | `1,101` | y-direction grid start/end indices `TAR_Y1, TAR_Y2` |
+| ③    | `1,50`  | z-direction grid start/end indices `TAR_Z1, TAR_Z2` |
+| ④    | `1e-5`  | Conductivity of this block `TAR_CONDUCTIVITY` (S/m) |
+
+> This example uses two "anomalous bodies" to assemble a half-space: block 1 z=1~50 (air, 1e-5 S/m) + block 2
+> z=51~100 (stratum, 1e-2 S/m). After the receiver count `Point_Num`, each measurement point occupies 2 lines
+> (index + coordinates relative to the source center).
+
+### 3.2 Terrain Mesh File Complex_Terrain
+
+The terrain is described by a **surface triangular
+mesh**, supported in two formats; keep **only one of them** in the folder;
+if both exist at the same time, the program takes `.dat` as the priority and warns that `.stl` is ignored.
+
+**Format 1: `Complex_Terrain.dat` (original text format)**
+
+```
+Number of Nodes and Elements:
+10039               ← total number of nodes n_point
+5426                ← total number of triangular elements n_face
+Nodes Coordinates:
+1  -21000.0  -21000.0  224.08   ← n_point lines: node index, X, Y, Z
+2  -21000.0  -20001.8  224.08
+...                        (a `!` comment may appear in a line)
+END Nodes Coordinates
+NormalAreaElements:
+1  1  2  10039       ← n_face lines: element index, node1, node2, node3
+...
+END NormalAreaElements
+```
+
+| Line                       | Content                                                                                                        |
+| -------------------------- | -------------------------------------------------------------------------------------------------------------- |
+| Line 1                     | Title line, arbitrary                                                                                          |
+| Line 2                     | Total number of nodes `n_point`                                                                                |
+| Line 3                     | Total number of triangular elements `n_face`                                                                   |
+| Line 4                     | Title line, arbitrary                                                                                          |
+| Lines 5 ~ 4+n_point        | One line per node: `node index, X, Y, Z`                                                                       |
+| The following 1 line       | Section end marker `END Nodes Coordinates` (skipped by the program as a title line)                            |
+| The following 1 line       | Face-section title `NormalAreaElements:` (skipped by the program as a title line)                              |
+| The following n_face lines | One line per element: `element index, node1 index, node2 index, node3 index` (nodes in counterclockwise order) |
+| The last 1 line            | End marker `END NormalAreaElements` (not read by the program)                                                  |
+
+> The two section-marker lines after the node section (GiD export) are in the same position as the "2 title lines" of the old format;
+> the program always skips them as title lines, so both notations are compatible.
+
+**Format 2: `Complex_Terrain.stl` (ASCII STL format)**
+
+The STL file format is a file format used to describe the surface geometry of three-dimensional objects, widely used in rapid prototyping, 3D printing, and computer-aided manufacturing (CAM).
+An STL file subdivides an object's surface into a series of small triangles, each defined by a normal vector and three vertex coordinates.
+STL files have two formats: text format (ASCII) and binary format.
+
+Standard ASCII STL uses the `facet`/`endfacet` element keywords,
+with nodes represented by `vertex` lines, for example in the following format:
+
+```
+facet normal nx ny nz
+    outer loop
+        vertex v1x v1y v1z
+        vertex v2x v2y v2z
+        vertex v3x v3y v3z
+    endloop
+endfacet
+```
+
+Users can create and edit STL files with common CAD software (such as AutoCAD, Blender, FreeCAD, MeshLab, SketchUp, Gid, Maya, 3ds Max, etc.).
+
+The program automatically handles two things when reading (no user action needed):
+
+1. **Vertex deduplication/merging**: in STL, each face writes its vertices independently; duplicate vertices (tolerance 1e-5) are automatically merged into a unique node table;
+2. **Orientation correction**: compare each face's cross-product direction with the `facet normal` in the file; if opposite, swap the face's 2nd and 3rd
+   nodes to guarantee that the normal-direction convention is consistent with the `.dat` format.
+
+### 3.3 Anomalous Body Mesh File Complex_anomalous
+
+Exactly the same as the terrain file: either the `.dat` or `.stl` format (when both exist, `.dat` takes priority),
+with the same reading, deduplication, and orientation-correction rules; the `.dat` section marker lines (`END Nodes Coordinates`,
+`NormalAreaElements:`, `END NormalAreaElements`) are the same as the terrain file and are likewise compatible with the program.
+
+The `Complex_anomalous.dat` / `.stl` in this directory describe **a complex three-dimensional anomalous body under undulating terrain**
+(2663 nodes, 5322 triangular elements; extent x ≈ -302 ~ 248 m, y ≈ -197 ~ 176 m,
+z ≈ 51 ~ 285 m, embedded near the terrain surface).
+Its conductivity is specified by line 27 of `input.dat`,
+`tao_abnormal = 4.0` (a low-resistivity body).
+
+> Mesh generation suggestion: use professional preprocessing software to build the terrain surface / anomalous-body surface triangular mesh and export it,
+> or choose "Export → STL" to generate an ASCII STL file.
+
+### 3.4 Modeling Notes
+
+1. **The number of grids in the x and y directions is recommended to be odd**, so that the model center (the source center) falls exactly at the center of a Yee grid face;
+   since the magnetic induction intensity B is defined at the grid-face center, the `dBz/dt` measurement point should **preferably be placed at the center of a Yee grid face**,
+   to guarantee the correspondence between the measurement-point response and the field definition.
+2. **The transmitting loop side length should be an odd multiple of the grid size**, so that the loop center falls on a Yee grid edge, ensuring strict alignment between the source-current edges and the grid edges.
+3. **When using the Dirichlet boundary (`Logic_PML=0`)**, the grid refinement region (core uniform area) should **cover the transmitting source,
+   the receiver points, and the anomalous body**, to guarantee the computational accuracy of the above region; the outer large grids are used to extend the computational domain and reduce
+   the influence of the zero-field boundary on the result.
+4. **The z-direction position of the transmitting source is by default at `NZS+1`** (`NZS=NZ/2`, i.e. the middle of the grid), determined automatically by the program,
+   and does not need to be specified in the input file.
+5. **Airborne (TEM) simulation**: the program also supports airborne scenarios — set the upper region to air, and
+   add several more layers of air grids below the source plane, so that both the transmitting source and the receiver points are in the air.
+
+---
+
+## 4. CPML Absorbing Boundary
+
+When `Logic_PML=1`, the program uses **CPML (Convolutional Perfectly Matched Layer, Roden & Gedney 2000)**
+as the absorbing boundary, absorbing the outward-propagating electromagnetic field at the periphery of the computational region to simulate an "infinite earth stratum"
+and avoid boundary reflections contaminating the late-time response.
+
+### 4.1 Implementation Location
+
+| File                                            | Function                                                                                                                                                                                                                                       |
+| ----------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `module/pml-parameters.f90`                     | Parameter declaration: σ/α/κ maxima, PML layer count, the 26 memory-variable ψ arrays, the convolution coefficients `b_e/c_e`, the scaling arrays `den_*(=1/κ)`, and the z-direction recursive coefficient `c_h_zz` for Hz                     |
+| `lib/get-pml-paramters.f90`                     | Builds the spatial distributions of σ/α/κ on the six x/y/z boundary faces (polynomial increasing from inside to outside, E/H staggered sampling), the `den_*` scaling arrays, and the convolution coefficients; called only when `Logic_PML=1` |
+| `lib/Iteration.f90` (subroutine Iteration_cpml) | In the field-update main loop, **embeds** the recurrence and correction terms of the memory variables ψ (24 ψ for Ex/Ey/Ez and Hx/Hy); the z direction of Hz uses a recursive convolution based on `c_h_zz` (not ψ)                            |
+
+### 4.2 Parameters and Their Meanings
+
+| Parameter           | Default value                               | Meaning                                                                                       |
+| ------------------- | ------------------------------------------- | --------------------------------------------------------------------------------------------- |
+| `PML_X,PML_Y,PML_Z` | line 5 of input.dat (≥5 layers recommended) | PML layer count in the three directions                                                       |
+| `ma`                | 3                                           | Polynomial order of σ along the thickness direction (power-law growth from inside to outside) |
+| `mb`                | 1                                           | Polynomial order of α along the thickness direction                                           |
+| `sig_max`           | 1.0e2                                       | Maximum conductivity at the outer side of the PML (determines absorption strength)            |
+| `alpha_max`         | **1.0e-1**                                  | Maximum value of the complex-frequency-shift factor                                           |
+| `kappa_max`         | 1.0                                         | Maximum value of the coordinate-stretching factor (1 means no stretching)                     |
+
+The spatial distributions of σ and α along the thickness (taking the lower side of the x direction as an example; the other boundaries are symmetric):
+
+```
+σ(i) = sig_max · ((L−i)/(L−1))^ma
+α(i) = alpha_max · ((i−1)/(L−1))^mb
+```
+
+The E field is sampled over whole layers and the H field over half layers (staggered), so the H-direction σ/α/κ are constructed with a half-layer offset
+(`i−0.5`), staggered with respect to the E direction.
+
+### 4.3 Usage Notes
+
+1. In `input.dat`, set line 4 `Logic_PML=1`, and give the PML layer counts for the three directions in line 5
+   (e.g. `15,15,15`); the grid size inside the PML layers should be consistent with the core area (kept uniform).
+2. **Taking α_max = 0.1 is the key tuning parameter**: α (the complex-frequency-shift factor) is responsible for absorbing the low-frequency diffusion field.
+   If it is too small (e.g. 0.01), the low-frequency reflected field in the late-time stage cannot decay in time and will bounce back and forth across the boundary,
+   causing exponential divergence on the order of 10⁻⁵ s after turn-off (the result becomes NaN) — this is
+   the most common source of instability in the CPML version, so be sure to keep α_max = 1.0e-1.
+3. To switch back to the original boundary: change line 4 to `0`; the behavior is exactly the same as the old version, no recompilation needed.
+
+### 4.4 Comparison with the Dirichlet Boundary
+
+|                    | CPML (`Logic_PML=1`)                                       | Dirichlet (`Logic_PML=0`)                                   |
+| ------------------ | ---------------------------------------------------------- | ----------------------------------------------------------- |
+| Boundary handling  | Absorbing layer, simulates an unbounded space              | Field is zero at the boundary                               |
+| Late-time accuracy | Absorbs reflections, decay curve is flat                   | Boundary reflections may contaminate the late-time response |
+| Computational cost | 24 extra memory-variable recurrences per step (about +30%) | No extra overhead                                           |
+| Grid requirement   | Uniform grid required inside the PML layers                | No special requirement                                      |
+| Stability          | Stable when parameters are tuned correctly                 | Stable                                                      |
+
+> **Consistency verification**: in a uniform half-space example (81×81×80 grid, 15 PML layers, 1 ms plateau,
+> compared after turn-off), the CPML and Dirichlet boundaries give consistent response curves in the early stage after turn-off (when the field has not yet reached the boundary),
+> with a difference < 0.2% (arising from floating-point rounding accumulation due to the different loop orderings in the two versions,
+> not a physical difference), proving that the CPML implementation is equivalent to and correct with respect to the main iteration.
+
+---
+
+## 5. Output File Description
+
+| File                            | Content                                                                                                                                                                                                                                         |
+| ------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| `dBzdt_1.txt`, `dBzdt_2.txt`, … | One file per receiver point. The first 2 lines of the header are descriptions (point index, point coordinates); each subsequent line has 3 columns: iteration step, time after turn-off (s), and the magnetic-field response value at that time |
+| `CTIME_TIXING_UPCOS.DAT`        | Computation time series; each line has 3 columns: cumulative time, time step, source current amplitude (waveform)                                                                                                                               |
+| `conductivity.vtk`              | Model conductivity distribution (regular-grid VTK format), can be opened with ParaView/Tecplot etc. to check that the model was built correctly                                                                                                 |
+| `logfile.log`                   | Run log: parameter echo, format-selection hints, run errors, etc.                                                                                                                                                                               |
+| `fort.5141`                     | Debug output of the conformal mesh computation process                                                                                                                                                                                          |
+| `TEM_decay_curve.png`           | Decay-curve figure (generated by `TEM_decay_plot.py` in [Section 6](#6-quick-decay-curve-plotting-tem_decay_plotpy))                                                                                                                            |
+
+At the end of the run the screen prints the iteration progress of each segment and the total computation time; on normal completion, the end of `logfile.log`
+shows `Computation finished!`.
+
+---
+
+## 6. Quick Decay-Curve Plotting (TEM_decay_plot.py)
+
+The `TEM_decay_plot.py` in the program directory is used to quickly plot the forward-modeling results `dBzdt_*.txt` as
+**decay voltage curves** (log-log coordinates).
+
+### 6.1 Usage
+
+```bash
+# requires numpy and matplotlib
+pip install numpy matplotlib
+
+# run in the directory where the computation output files (dBzdt_*.txt) are located
+python TEM_decay_plot.py
+```
+
+The script automatically searches for all `dBzdt_*.txt` files in the same directory, draws one curve per receiver point,
+and automatically annotates the point index and coordinates in the legend; by default it outputs a high-resolution image **`TEM_decay_curve.png`**
+(dpi=600) and pops it up.
+
+### 6.2 Main Adjustable Parameters (the "User parameters" section at the top of the script)
+
+| Parameter              | Default value         | Description                                                                                                               |
+| ---------------------- | --------------------- | ------------------------------------------------------------------------------------------------------------------------- |
+| `file_pattern`         | `dBzdt_*.txt`         | Pattern for matching result files                                                                                         |
+| `xmin, xmax`           | `1e-6, 1e-1`          | Display range of the horizontal axis (time, s)                                                                            |
+| `ymin, ymax`           | `None, None`          | Display range of the vertical axis (response); `None` means automatic                                                     |
+| `use_abs`              | `True`                | `True` plots \|dBz/dt\| (positive response), `False` plots the signed value; the vertical-axis label switches accordingly |
+| `savefig`              | `True`                | Whether to save the figure                                                                                                |
+| `save_name`            | `TEM_decay_curve.png` | Saved file name                                                                                                           |
+| `dpi`                  | `600`                 | Image resolution                                                                                                          |
+| `label_fontsize`, etc. | 18/15/15              | Font sizes for labels, ticks, and legend                                                                                  |
+| `linewidth`            | `2.5`                 | Curve line width                                                                                                          |
+
+When there are more than 8 curves, a sequential blue-gradient colormap is automatically used, so colors are not recycled.
+
+> Curve value rule: only points with `time > 0` and response value > 0 are plotted (negative/zero values cannot be displayed in log-log coordinates).
+
+---
+
+## 7. Quick Start (the example bundled in this directory)
+
+This directory bundles an example: **a complex three-dimensional anomalous body under undulating terrain** (model 101×101×100 grids, grid size
+20 m, source side length 500 m, background 0.01 S/m, the anomalous body is a low-resistivity body of 4.0 S/m, 1 measurement point (directly below
+the source center); anomalous-body mesh 2663 nodes/5322 elements, terrain mesh 10039 nodes/5426 elements).
+
+Run steps:
+
+1. Confirm that the directory contains: `input.dat`, `Complex_anomalous.dat` (and `.stl`, describing the same anomalous body),
+   `Complex_Terrain.dat` (and `.stl`, describing the same terrain). When the `.dat` and `.stl` of the same kind
+   both exist, the program takes `.dat` as the priority and gives a hint.
+2. Open `tem3dfdtd.sln` in VS2019 → select `Release | x64` → build.
+3. Copy the generated exe to this directory (or put the input files in the exe directory) and run.
+4. Observe the screen output; the normal flow is:
+   ```
+   Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.
+   Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.
+   The number of grids in the core area is odd
+   At least 320M memory is needed!
+   ...
+   Conformal mesh of terrain is complete!
+   Ray tracing computation of terrain is complete!
+   Conformal mesh of terrain is finished
+   Start conformal processing of the anomalous body
+   ...
+   Now computing fraction:  1
+   50 steps have just finished
+   ...
+   ```
+5. After the computation finishes, check the output files `dBzdt_1.txt` and `conductivity.vtk`.
+6. (Optional) Run `python TEM_decay_plot.py` to generate the decay-curve figure `TEM_decay_curve.png`
+   (see [Section 6](#6-quick-decay-curve-plotting-tem_decay_plotpy)).
+
+---
+
+## 8. FAQ
+
+**Q1: Running reports `libiomp5md.dll` not found**
+The OpenMP runtime library is missing. Copy Intel oneAPI's `bin/libiomp5md.dll` next to the exe
+(or add it to PATH).
+
+**Q2: The prompt `Both Complex_Terrain.dat and Complex_Terrain.stl exist! ...` appears**
+Both format files are present. The program gives `.dat` priority. If you want to use STL, move the `.dat` file away or rename it.
+
+**Q3: The computation is very slow / out of memory**
+Reduce `NX,NY,NZ` or increase `GridSize`; control `NSTOP`; `MAX_OFF_TIME` determines the actual number of iteration steps,
+and the program uses the smaller of the two. Before running, the required memory estimate is printed.
+
+**Q4: How to compute only a uniform half-space (no anomalous body, no terrain)?**
+A half-space model should contain two parts: "air + earth". Move `Complex_anomalous.*` and
+`Complex_Terrain.*` away, and set `TEMP_II = 2` in `input.dat`: set block 1 to the upper part
+(air, conductivity e.g. `1e-5`), block 2 to the lower part (earth, conductivity e.g. `0.01`), which forms a
+uniform half-space. Note: **when `TEMP_II = 0`, the whole model is filled only with the background conductivity (a whole-space uniform medium,
+without an air layer)**.
+
+**Q5: How to write the receiver coordinates?**
+The coordinates are **local coordinates relative to the loop-source center** (unit: m), with positive/negative directions consistent with the coordinate axes.
+
+**Q6: How to choose between the CPML absorbing boundary and the original Dirichlet boundary?**
+Line 4 of `input.dat` is the switch `Logic_PML`: `1` enables the CPML absorbing boundary (line 5
+`10,10,10` is the PML layer count in the three directions, which you can adjust), which effectively absorbs boundary reflections and gives a flatter late-time
+(large-offset / late-time) decay curve; `0` uses the original non-uniform-grid Dirichlet (zero-field)
+boundary. Switching does not require recompilation.
+
+---
+
+## 9. References
+
+[1] Sun Huaifeng, Li Xiu, Li Shucai, et al. Three-dimensional FDTD forward modeling of loop-source TEM considering the turn-off time [J]. Chinese Journal of Geophysics, 2013, 56(3): 1049-1064.
+
+[2] Liu Shangbin, Li Xuefeng, Lan Riyan, et al. CPML absorbing boundary for the low-frequency-approximated Maxwell equations in transient electromagnetics and its application method [J]. Chinese Journal of Geophysics, 2022, 65(4): 1472-1481.
+
+[3] Li X, Zhao Q, Hu S, et al. Introducing complex geometries to Yee cells in FDTD for transient electromagnetic forward modeling[J]. Geophysics, 2025, 91(2): F1-F12.
+
+## 10. Contributors
+
+The entire project was carried out under the leadership of Professor Sun Huaifeng of Shandong University. The main contributors are already noted in the code comments and references. To contact us, please visit https://faculty.sdu.edu.cn/sun/
+
+In addition to this open-source code repository, we also provide a commercial-version software or dedicated solver that supports efficient GPU computation. If needed, please visit https://em3d.cn for more information.
+
+## 11. Statement
+
+- For contact, please use sunhuaifeng@email.sdu.edu.cn, and do not continue using the gmail address noted in the code, because that gmail address often has sending/receiving problems. Thank you.
+- The tdem.org listed in the code cannot be maintained for the time being due to limited energy.
+- All code is maintained and released through the git repository https://git.em3d.cn/
@@ -1,9 +1,7 @@
-# tem3dfdtd
-
-请访问https://em3d.cn/ 获取更多信息。
-
 # 三维瞬变电磁 FDTD 正演程序 (tem3dfdtd)

+[中文](./README.md) | [English](./README.en.md)
+
 本程序基于 **FDTD（时域有限差分）方法**对三维瞬变电磁（TEM）响应进行正演模拟。核心算法采用 Wang–Hohmann (1993)改进的 Du Fort–Frankel 方法：在 Yee 网格上对磁场 H 迭代求解，通过引入**虚介电常数**保证显式迭代的时间稳定性；采用**共形网格技术**处理起伏地形与任意形状异常体（以表面三角网格描述）；支持矩形回线源，可计算地面 TEM 、半航空（SATEM）全航空（ATEM）等模式。

 程序的整体框架与三大核心技术分别源于以下工作：
@@ -20,7 +18,6 @@

 ## 目录

- [tem3dfdtd](#tem3dfdtd)
 - [三维瞬变电磁 FDTD 正演程序 (tem3dfdtd)](#三维瞬变电磁-fdtd-正演程序-tem3dfdtd)
  - [目录](#目录)
  - [1. 运行环境与编译](#1-运行环境与编译)
@@ -79,7 +76,7 @@
 3. 选择 **`Release | x64`** 配置,点击 **生成 → 生成解决方案**。
 4. 运行前准备:程序在工作目录中查找 `input.dat` 及网格文件,因此请将
   `input.dat`、`Complex_Terrain.dat/.stl`、`Complex_anomalous.dat/.stl` 放到
-   **exe 所在目录**(或通过"项目属性 → 调试 → 工作目录"指定)。
+   **tem3dfdtd-open\tem3dfdtd**目录中(或通过"项目属性 → 调试 → 工作目录"指定)。
 5. 直接运行 `tem3dfdtd.exe`,或按 F5 调试运行。

 > **注意**:程序使用 OpenMP,运行时需要 Intel 的运行时库 `libiomp5md.dll`
@@ -538,7 +535,6 @@ OpenMP 运行时库缺失。将 Intel oneAPI 安装目录 `bin/libiomp5md.dll`

 ## 11.声明

- 如需联系孙怀凤教授，请使用sunhuaifeng@email.sdu.edu.cn，不要继续使用代码中标注的gmail邮箱了，因为gmail邮箱经常会出现收发问题，谢谢。
- 代码中列出的sunhuaifeng.com网站域名已经注销，不再属于课题组，也已经永久不再维护，后期将会修正代码，请勿访问
+- 如需联系，请使用sunhuaifeng@email.sdu.edu.cn，不要继续使用代码中标注的gmail邮箱了，因为gmail邮箱经常会出现收发问题，谢谢。
 - 代码中列出的tdem.org 由于精力原因暂时无法维护
 - 所有代码均通过git仓库进行维护和发布https://git.em3d.cn/
@@ -1,4 +1,4 @@
-
+
 Microsoft Visual Studio Solution File, Format Version 12.00
 # Visual Studio Version 16
 VisualStudioVersion = 16.0.35826.135
@@ -1,552 +0,0 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-
-# 三维瞬变电磁 FDTD 正演程序 (tem3dfdtd-cmm)
-
-本程序基于 **FDTD（时域有限差分）方法**对三维瞬变电磁（TEM）响应进行正演模拟。核心算法采用 Wang–Hohmann (1993)改进的 Du Fort–Frankel 方法：在 Yee 网格上对磁场 H 迭代求解，通过引入**虚介电常数**保证显式迭代的时间稳定性；采用**共形网格技术**处理起伏地形与任意形状异常体（以表面三角网格描述）；支持矩形回线源，可计算地面 TEM 、半航空（SATEM）全航空（ATEM）等模式。
-
-程序的整体框架与三大核心技术分别源于以下工作：
-
- **整体框架与核心迭代算法**（考虑关断时间的回线源激发 TEM 三维时域有限差分正演，Wang–Hohmann 改进的 Du Fort–Frankel 方法、虚介电常数、含关断时间的源波形）[1]；
- **CPML 吸收边界**（瞬变电磁低频近似 Maxwell 方程的 CPML 吸收边界及施加方法）[2]；
- **共形网格技术**（Introducing complex geometries to Yee cells in FDTD for transient electromagnetic forward modeling）[3]。
-
-各方法的具体原理、公式推导与实现细节详见文献 [1]–[3]，本文档不做方法层面的展开介绍。
-
-代码结构：`main.f90`（主程序）、`module/`（全局参数与模块）、`lib/`（各功能子程序）。
-
---
-
-## 目录
-
- [三维瞬变电磁 FDTD 正演程序 (tem3dfdtd-cmm)](#三维瞬变电磁-fdtd-正演程序-tem3dfdtd-cmm)
-  - [目录](#目录)
-  - [1. 运行环境与编译](#1-运行环境与编译)
-    - [1.1 环境要求](#11-环境要求)
-    - [1.2 用 VS2019 编译运行(推荐)](#12-用-vs2019-编译运行推荐)
-    - [1.3 命令行编译(可选)](#13-命令行编译可选)
-    - [1.4 Linux 环境编译(makefile)](#14-linux-环境编译makefile)
-  - [2. 程序流程与模块结构](#2-程序流程与模块结构)
-    - [模块 1:程序控制与参数输入](#模块-1程序控制与参数输入)
-    - [模块 2:Yee 网格生成](#模块-2yee-网格生成)
-    - [模块 3:电性参数构建](#模块-3电性参数构建)
-    - [模块 4:激励源与时间序列](#模块-4激励源与时间序列)
-    - [模块 5:FDTD 电磁场计算](#模块-5fdtd-电磁场计算)
-    - [模块 6:输出](#模块-6输出)
-  - [3. 输入文件格式](#3-输入文件格式)
-    - [3.1 参数控制文件 input.dat](#31-参数控制文件-inputdat)
-    - [3.2 地形网格文件 Complex_Terrain](#32-地形网格文件-complex_terrain)
-    - [3.3 异常体网格文件 Complex_anomalous](#33-异常体网格文件-complex_anomalous)
-    - [3.4 建模注意事项](#34-建模注意事项)
-  - [4. CPML 吸收边界](#4-cpml-吸收边界)
-    - [4.1 实现位置](#41-实现位置)
-    - [4.2 参数与含义](#42-参数与含义)
-    - [4.3 使用说明](#43-使用说明)
-    - [4.4 与 Dirichlet 边界的对比](#44-与-dirichlet-边界的对比)
-  - [5. 输出文件说明](#5-输出文件说明)
-  - [6. 衰减曲线快速成图(TEM_decay_plot.py)](#6-衰减曲线快速成图tem_decay_plotpy)
-    - [6.1 使用方式](#61-使用方式)
-    - [6.2 主要可调参数(脚本顶部"User parameters"区)](#62-主要可调参数脚本顶部user-parameters区)
-  - [7. 快速上手(本目录自带算例)](#7-快速上手本目录自带算例)
-  - [8. 常见问题](#8-常见问题)
-  - [9. 参考文献](#9-参考文献)
-
---
-
-## 1. 运行环境与编译
-
-### 1.1 环境要求
-
-| 项目           | 要求                                                                                                             |
-| -------------- | ---------------------------------------------------------------------------------------------------------------- |
-| 操作系统       | Windows 10/11(64 位)或 Linux(x86_64)                                                                             |
-| 集成环境       | Windows:Visual Studio 2019;Linux:make                                                                            |
-| Fortran 编译器 | Windows:Intel oneAPI Fortran(ifort/ifx); Linux:NVIDIA HPC SDK(nvfortran),见 [1.4 节](#14-linux-环境编译makefile) |
-| 并行支持       | OpenMP(多核 CPU 加速,代码中已使用 `!$OMP PARALLEL DO` 指令)                                                      |
-
-### 1.2 用 VS2019 编译运行(推荐)
-
-1. 安装 Visual Studio 2019(勾选"C++ 桌面开发"工作负载)与 Intel oneAPI(安装时勾选
-   "Intel Fortran Compiler 的 Visual Studio 集成")。
-2. 双击打开工程文件 **`tem3dfdtd-cmm.sln`**(Intel Fortran 工程,对应
-   `tem3dfdtd-cmm.vfproj`)。工程已包含全部源文件,配置说明:
-   - `Debug | x64` / `Release | x64`:使用 **ifx** 编译器(oneAPI 默认),**推荐**;
-   - `Debug | Win32` / `Release | Win32`:使用 ifort 编译器。
-3. 选择 **`Release | x64`** 配置,点击 **生成 → 生成解决方案**。
-4. 运行前准备:程序在工作目录中查找 `input.dat` 及网格文件,因此请将
-   `input.dat`、`Complex_Terrain.dat/.stl`、`Complex_anomalous.dat/.stl` 放到
-   **exe 所在目录**(或通过"项目属性 → 调试 → 工作目录"指定)。
-5. 直接运行 `tem3dfdtd-cmm.exe`,或按 F5 调试运行。
-
-> **注意**:程序使用 OpenMP,运行时需要 Intel 的运行时库 `libiomp5md.dll`
-> (位于 Intel oneAPI 安装目录 `bin/` 下)。若提示缺少该 DLL,可将它复制到
-> exe 同目录(本目录已放置一份),或将其所在目录加入系统 PATH。
-
-### 1.3 命令行编译(可选)
-
-在"Intel oneAPI 命令行"环境(oneAPI Command Prompt)下:
-
-```bat
-ifx -c -O2 -Qopenmp module\*.f90
-ifx -c -O2 -Qopenmp -Qopenmp lib\*.f90 main.f90
-ifx -O2 -Qopenmp *.obj -o tem3dfdtd-cmm.exe
-```
-
-(链接时需要 MSVC 的 link.exe 与 Windows SDK 库,建议直接使用 VS 的"开发人员命令提示符",
-并在 PATH 中加入 Intel oneAPI 的 `bin` 目录。)
-
-### 1.4 Linux 环境编译(makefile)
-
-程序根目录提供 `makefile`,用于在 Linux 系统上编译,使用的编译器为
-**NVIDIA HPC SDK**(命令名 `nvfortran`)。makefile 内容如下:
-
-```makefile
-obj:
-	nvfortran -fast -g -traceback   -mp=multicore -Minfo=mp -c module/*.f90 lib/*.f90 main.f90
-exe:
-	nvfortran -fast -g -traceback   -mp=multicore -Minfo=mp -o main.exe *.o
-clean:
-	 rm -rf *.exe *.mod *.o nohup.out *.txt *.log SUR*  *.DAT *.obj
-```
-
-**编译步骤**
-
-```bash
-# 1. 加载 NVIDIA HPC SDK(以 module 方式安装时)
-module load nvhpc        # 或自行将 nvhpc 的 bin/ 加入 PATH
-
-# 2. 编译(两个目标按顺序执行;make 默认只执行第一个目标 obj,故必须再执行 exe)
-make obj && make exe
-
-# 3. 运行(先把 input.dat 及网格文件放到 exe 所在目录)
-./main.exe
-```
-
-**说明与注意事项**
-
-1. **目标结构**:makefile 的 `obj`/`exe` 两个目标之间没有依赖关系,`make` 单独执行
-   只会编译出 `.o` 文件,**必须依次执行 `make obj && make exe`** 才能得到 `main.exe`。
-2. **OpenMP 线程数**:默认使用全部逻辑核,可运行前用 `export OMP_NUM_THREADS=N` 限制。
-3. **`make clean` 会一并删除计算结果**(`*.txt`、`*.log`、`*.DAT`、`nohup.out`、
-   `SUR*` 等),执行前请先备份 `dBzdt_*.txt` 等输出文件。
-4. 与 Windows 版的差异:Linux 下可执行文件名为 `main.exe`,工作目录中同样需要
-   `input.dat` 与 `Complex_Terrain.*`、`Complex_anomalous.*` 网格文件。
-
---
-
-## 2. 程序流程与模块结构
-
-主程序执行顺序(`main.f90`):
-
-```
-GETDATA → CHECKPARAMETERS → MEMORY_USE_ESTIMATION → ALLOCATEMEMORY
-→ GET_NON_UNIFORMGRID → ZERO → GET_COORDINATES
-→ Get_Receiver_Gridlabel → RES_CONFIGURE → TIME_SERIOUS
-→ Get_eps_r → (Logic_PML=1 时) Get_pml_parameters → Get_mstop
-→ GetSourcePosition → Iteration → FREE_MEMORY
-```
-
-按功能划分为以下 6 个模块:
-
-### 模块 1:程序控制与参数输入
-
-| 文件                            | 功能                                                                       |
-| ------------------------------- | -------------------------------------------------------------------------- |
-| `main.f90`                      | 主程序,控制整个计算流程                                                    |
-| `lib/getdata.f90`               | 读取参数控制文件 `input.dat`;检测地形/异常体网格文件的存在性并选择读取格式 |
-| `lib/checkparameters.f90`       | 将读入的计算参数回显到 `logfile.log`,便于人工检查                          |
-| `lib/memory-use-estimation.f90` | 根据网格规模估算所需内存并打印提示                                         |
-
-### 模块 2:Yee 网格生成
-
-| 文件                          | 功能                                                                       |
-| ----------------------------- | -------------------------------------------------------------------------- |
-| `lib/allocatememory.f90`      | 根据输入参数动态分配所有全局数组(含 CPML 记忆变量数组,仅 `Logic_PML=1` 时) |
-| `lib/get_non_uniformgrid.f90` | 生成 x/y/z 三方向的非均匀网格(核心区均匀 + 外围按 1.3 倍递增扩展)          |
-| `lib/get_coordinates.f90`     | 计算各网格节点(含 Yee 节点)的坐标,坐标以源中心为原点                       |
-| `lib/zero.f90`                | 将所有电磁场数组初始化为 0;den\_\* 置 1、c_h_zz 置 0、CPML 记忆变量清零    |
-
-### 模块 3:电性参数构建
-
-| 文件                                | 功能                                                                                                                                                                           |
-| ----------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
-| `lib/resistivity-configuration.f90` | 构建模型电导率:无地形时按背景电导率 + 块状异常体赋值;有地形时调用共形网格;最后将电导率分配到 x/y/z 三个方向的棱边并写出 `conductivity.vtk`                                     |
-| `lib/Terrain_conformal.f90`         | **地形共形网格**:从 `Complex_Terrain.dat/.stl` 读入地形三角网格,用射线–三角形求交(Möller–Trumbore 算法)沿 x/y/z 三个方向填充每个棱边的等效电导率,处理起伏地形与空气/地层的分界 |
-| `lib/Anomalous_conformal.f90`       | **异常体共形网格**:从 `Complex_anomalous.dat/.stl` 读入异常体表面网格,采用同样的射线求交方法将异常体电导率(`tao_abnormal`)填充到棱边                                           |
-
-> **说明**:存在地形文件时,地形分支中空气电导率取 `AIR_CONDUCTIVITY = 1e-6 S/m`,
-> 地层电导率取**块状异常体参数中的第 2 个电导率** `TAR_CONDUCTIVITY(2)`(见 input.dat
-> 第 9 行起的第二组数据);此时 input.dat 中块状异常体本身不直接生效,而是以 `Complex_anomalous`
-> 网格文件描述异常体、以 `tao_abnormal` 赋予其电导率。
-
-### 模块 4:激励源与时间序列
-
-| 文件                          | 功能                                                                                                      |
-| ----------------------------- | --------------------------------------------------------------------------------------------------------- |
-| `lib/time-serious.f90`        | 生成整个计算的时间序列(含源波形),并根据 `MAX_OFF_TIME` 校正迭代步数 `NSTOP`;写出 `CTIME_TIXING_UPCOS.DAT` |
-| `lib/tixing-source-upcos.f90` | 源波形:**梯形 + 余弦上升**的关断电流波形(常用,`SOURCE_TYPE = 'TIXING_UPCOS'`)                             |
-| `lib/tixing-source.f90`       | 纯梯形波形源(`TIXING_RAMP`)                                                                               |
-| `lib/sin-source.f90`          | 半正弦波形源(`HALF_SIN`)                                                                                  |
-| `lib/triangle-source.f90`     | 三角波形源(`TRIANGLE`)                                                                                    |
-| `lib/get-eps-r.f90`           | 计算虚介电常数 `EPS_R = 3·(Δt/Δx)²/μ₀` 及迭代系数,保证显式 FDTD 稳定                                      |
-| `lib/get-mstop.f90`           | 将总迭代过程切成若干"计算分段",每段独立分配缓存,便于内存管理                                              |
-
-### 模块 5:FDTD 电磁场计算
-
-| 文件                        | 功能                                                                                                                                                                                                                                                                                  |
-| --------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
-| `lib/GetSourcePosition.f90` | 确定回线源在网格中的位置,标记源所在的棱边(电流赋值区域)                                                                                                                                                                                                                               |
-| `lib/Iteration.f90`         | **核心迭代子程序**:按分段循环推进时间步,更新 Ex/Ey/Ez、Hz 场,按源波形加载电流;在每个分段末尾对各接收点计算 Hz(由环绕该点的 8 个节点 Ex、Ey 差商加权得到)并写入结果文件。边界条件按 `Logic_PML` 开关切换:`1` 时在每个场更新中附加 CPML 记忆变量修正,`0` 时恢复原始 Dirichlet(零场)边界 |
-| `module/pml-parameters.f90` | **CPML 吸收边界模块**(Roden–Gedney 卷积 PML):σ/α/κ 多项式缩放参数(ma=3、mb=1)、26 个记忆变量 ψ 数组、b*e/c_e 卷积系数数组与 den*\*(=1/κ)缩放数组的声明                                                                                                                                |
-| `lib/get-pml-paramters.f90` | 构建 x/y/z 六个边界面的 σ/α/κ 分布(多项式从边界向内衰减,E/H 交错采样)及 den\_\* 缩放数组;仅在 `Logic_PML=1` 时由 main 调用                                                                                                                                                            |
-
-### 模块 6:输出
-
-| 文件                                          | 功能                                                |
-| --------------------------------------------- | --------------------------------------------------- |
-| `lib/Iteration.f90`(输出部分)                 | 写出各接收点响应文件 `dBzdt_1.txt`、`dBzdt_2.txt` … |
-| `lib/resistivity-configuration.f90`(输出部分) | 写出模型电导率分布 `conductivity.vtk`               |
-| `lib/time-serious.f90`(输出部分)              | 写出时间序列 `CTIME_TIXING_UPCOS.DAT`               |
-| `lib/free-memory.f90`                         | 计算结束后释放所有动态内存                          |
-
---
-
-## 3. 输入文件格式
-
-程序运行需要以下文件(全部放在 exe 的工作目录中):
-
-| 文件                              | 是否必须       | 说明                                      |
-| --------------------------------- | -------------- | ----------------------------------------- |
-| `input.dat`                       | **必须**       | 计算参数控制文件                          |
-| `Complex_anomalous.dat` 或 `.stl` | 可选(任选其一) | 异常体表面三角网格;缺省时模型视为均匀背景 |
-| `Complex_Terrain.dat` 或 `.stl`   | 可选(任选其一) | 地形表面三角网格;缺省时不考虑起伏地形     |
-
-### 3.1 参数控制文件 input.dat
-
-**自由格式读取**,按行顺序读取;数值后可加 `!` 注释(可整行注释或行尾注释)。
-下面以本目录自带的 `input.dat` 为例逐行说明:
-
-| 行号  | 示例           | 含义                                                                                      |
-| ----- | -------------- | ----------------------------------------------------------------------------------------- |
-| 1     | `1`            | 计算模式 `CAL_TYPE`:`1` = 地面 TEM,`2` = 半航空(SATEM)                                    |
-| 2     | `500`          | 发射回线边长 `SourceLength`(m)                                                            |
-| 3     | `101,101,100`  | x、y、z 三个方向的网格数 `NX,NY,NZ`                                                       |
-| 4     | `1`            | **边界条件开关 `Logic_PML`**:`1` = CPML 吸收边界,`0` = 原始非均匀网格 Dirichlet(零场)边界 |
-| 5     | `10,10,10`     | **PML 层数 `PML_X,PML_Y,PML_Z`**(x、y、z 方向,仅开关=1 时有效;建议 ≥ 5 层)                |
-| 6     | `25,25`        | x 方向:核心均匀网格起始/结束区间编号 `UniGridNumX1,UniGridNumX2`                          |
-| 7     | `25,25`        | y 方向:核心均匀网格区间编号 `UniGridNumY1,UniGridNumY2`                                   |
-| 8     | `20,30`        | z 方向:核心均匀网格区间编号 `UniGridNumZ1,UniGridNumZ2`                                   |
-| 9     | `20`           | 核心区均匀网格尺寸 `GridSize`(m)                                                          |
-| 10    | `0.01`         | 背景介质电导率 `BACKGROUND_CONDUCTIVITY`(S/m)                                             |
-| 11    | `2`            | 块状异常体数量 `TEMP_II`(无地形时按棱柱体填充;设 0 表示均匀模型)                          |
-| 12–15 | 见下           | **第 1 个块状异常体参数**,共 4 行                                                         |
-| 16–19 | 见下           | **第 2 个块状异常体参数**,共 4 行                                                         |
-| 20    | `4000000`      | 最大迭代次数 `NSTOP`                                                                      |
-| 21    | `90.101`       | 最大计算时间 `MAX_OFF_TIME`(**单位 ms**)                                                  |
-| 22    | `1e-6,1e-9`    | 上升沿持续时间与时间步 `RAISETIME, RAISESTEP`(s)                                          |
-| 23    | `60000e-6`     | 平台阶段持续时间 `WAVE`(s,即 60 ms)                                                       |
-| 24    | `1e-7,1e-9`    | 下降沿持续时间与时间步 `RAMP, RAMPSTEP`(s)                                                |
-| 25    | `1e-9`         | 初始时间步 `TIMESTEP`(s)                                                                  |
-| 26    | `1`            | 发射电流幅度 `AMP`(A)                                                                     |
-| 27    | `4.0`          | 异常体电导率 `tao_abnormal`(S/m,配合 `Complex_anomalous` 文件使用)                        |
-| 28    | `TIXING_UPCOS` | 源类型 `SOURCE_TYPE`:`TIXING_UPCOS` / `TIXING_RAMP` / `HALF_SIN` / `TRIANGLE`             |
-| 29    | `1`            | 接收点数 `Point_Num`                                                                      |
-| 30    | `1`            | 第 1 个接收点的编号                                                                       |
-| 31    | `0,0,0`        | 第 1 个接收点坐标(相对源中心,m)                                                           |
-| 32–33 | `2` / `x,y,z`  | (如有多余测点)第 2 个接收点(编号 + 坐标)                                                  |
-
-每个**块状异常体**由连续的 4 行组成:
-
-| 行  | 示例    | 含义                                |
-| --- | ------- | ----------------------------------- |
-| ①   | `1,101` | x 方向网格起止编号 `TAR_X1, TAR_X2` |
-| ②   | `1,101` | y 方向网格起止编号 `TAR_Y1, TAR_Y2` |
-| ③   | `1,50`  | z 方向网格起止编号 `TAR_Z1, TAR_Z2` |
-| ④   | `1e-5`  | 该块电导率 `TAR_CONDUCTIVITY`(S/m)  |
-
-> 本算例用两块"异常体"拼出半空间:第 1 块 z=1~50(空气,1e-5 S/m)+ 第 2 块
-> z=51~100(地层,1e-2 S/m)。接收点个数 `Point_Num` 后按每个测点 2 行排列
-> (编号 + 相对源中心坐标)。
-
-### 3.2 地形网格文件 Complex_Terrain
-
-地形由**表面三角网格**描述,支持两种格式,文件夹中**只保留其中一个**;
-若两个同时存在,程序以 `.dat` 为优先并提示 `.stl` 被忽略。
-
-**格式 1:`Complex_Terrain.dat`(原始文本格式)**
-
-```
-Number of Nodes and Elements:
-10039               ← 节点总数 n_point
-5426                ← 三角形单元总数 n_face
-Nodes Coordinates:
-1  -21000.0  -21000.0  224.08   ← n_point 行:节点编号, X, Y, Z
-2  -21000.0  -20001.8  224.08
-...                        (行中可带 ! 注释)
-END Nodes Coordinates
-NormalAreaElements:
-1  1  2  10039       ← n_face 行:单元编号, 节点1, 节点2, 节点3
-...
-END NormalAreaElements
-```
-
-| 行                  | 内容                                                                       |
-| ------------------- | -------------------------------------------------------------------------- |
-| 第 1 行             | 标题行,可任意                                                              |
-| 第 2 行             | 节点总数 `n_point`                                                         |
-| 第 3 行             | 三角形单元总数 `n_face`                                                    |
-| 第 4 行             | 标题行,可任意                                                              |
-| 第 5 ~ 4+n_point 行 | 每个节点一行:`节点编号, X, Y, Z`                                           |
-| 其后 1 行           | 区段结束标记 `END Nodes Coordinates`(程序按标题行跳过)                     |
-| 其后 1 行           | 面区标题 `NormalAreaElements:`(程序按标题行跳过)                           |
-| 其后 n_face 行      | 每个单元一行:`单元编号, 节点1编号, 节点2编号, 节点3编号`(节点按逆时针绕向) |
-| 末尾 1 行           | 结束标记 `END NormalAreaElements`(程序不读取)                              |
-
-> 节点区后的两个区段标记行(GiD 导出)与老版本格式的"2 行标题行"位置一致,
-> 程序一律按标题行跳过,因此两种写法均兼容。
-
-**格式 2:`Complex_Terrain.stl`(ASCII STL 格式)**
-
-标准 ASCII STL(如由 GiD、Gmesh 等导出),单元关键字为 `facet`/`endfacet`,
-节点用 `vertex` 行表示:
-
-```
-solid GiD
-  facet normal 0.000000e+00 0.000000e+00 1.000000e+00
-    outer loop
-      vertex  4.765000e+01 -1.318000e+01  9.550000e+00
-      vertex -1.318000e+01 -2.674000e+01 -1.544000e+01
-      vertex -1.003000e+01  4.387000e+01  2.476000e+01
-    endloop
-  endfacet
-  ...
-endsolid GiD
-```
-
-程序读取时自动处理两点(无需用户操作):
-
-1. **顶点去重合并**:STL 中每个面独立写顶点,重复顶点(容差 1e-5)自动合并为唯一节点表;
-2. **方向校正**:比较每个面的叉积方向与文件中的 `facet normal`,若相反则交换该面第 2、3
-   个节点,保证法线方向约定与 `.dat` 格式一致。
-
-### 3.3 异常体网格文件 Complex_anomalous
-
-与地形文件完全相同:`.dat` / `.stl` 两种格式任选其一(同时存在时 `.dat` 优先),
-读取方式、去重与方向校正规则均一致;`.dat` 的区段标记行(`END Nodes Coordinates`、
-`NormalAreaElements:`、`END NormalAreaElements`)与地形文件一致,同样与程序兼容。
-
-本目录的 `Complex_anomalous.dat` / `.stl` 描述的是**起伏地形下的一个复杂三维异常体**
-(2663 个节点、5322 个三角形单元;范围 x ≈ -302 ~ 248 m、y ≈ -197 ~ 176 m、
-z ≈ 51 ~ 285 m,嵌入地形面附近)。其电导率由 `input.dat` 第 27 行
-`tao_abnormal = 4.0` 指定(低阻体)。
-
-> 网格生成建议:用专业前处理软件(如 GiD)建立地形面/异常体表面三角形网格后导出,
-> 或选择"导出 → STL"生成 ASCII STL 文件。
-
-### 3.4 建模注意事项
-
-1. **x、y 方向网格数建议设为奇数**,使模型中心(源中心)恰好落在 Yee 网格面中心;
-   由于磁感应强度 B 定义于网格面中心,`dBz/dt` 测点应**优先布置在 Yee 网格面中心位置**,
-   以保证测点响应与场定义的对应关系。
-2. **发射回线边长应为网格尺寸的奇数倍**,使回线中心落在 Yee 网格棱边位置,保证
-   源电流棱边与网格棱边严格对位。
-3. **采用 Dirichlet 边界(`Logic_PML=0`)时**,网格加密区域(核心均匀区)应**覆盖发射源、
-   接收测点与异常体范围**,保证上述区域的计算精度,外围大网格用于扩展计算域、减弱
-   零场边界对结果的影响。
-4. **发射源的 z 方向位置默认在 `NZS+1` 处**(`NZS=NZ/2`,即网格中部),由程序自动
-   确定,无需在输入文件中指定。
-5. **全航空(TEM)模拟**:本程序同样支持全航空场景——将上半区设置为空气,并在
-   源平面以下再多布置若干层空气网格,即可保证发射源与接收测点均处于空气中。
-
---
-
-## 4. CPML 吸收边界
-
-本程序在 `Logic_PML=1` 时采用 **CPML(卷积完美匹配层, Roden & Gedney 2000)**
-作为吸收边界,在计算区域外围吸收向外传播的电磁场,模拟"无限大地层",
-避免边界反射污染晚时响应。
-
-### 4.1 实现位置
-
-| 文件                                       | 作用                                                                                                                            |
-| ------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------- |
-| `module/pml-parameters.f90`                | 参数声明:σ/α/κ 最大值、PML 层数、26 个记忆变量 ψ 数组、卷积系数 `b_e/c_e`、缩放数组 `den_*(=1/κ)`、Hz 的 z 向递归系数 `c_h_zz`  |
-| `lib/get-pml-paramters.f90`                | 构建 x/y/z 六个边界面内 σ/α/κ 的空间分布(多项式由内向外递增,E/H 交错采样)、`den_*` 缩放数组与卷积系数;仅在 `Logic_PML=1` 时调用 |
-| `lib/Iteration.f90`(子程序 Iteration_cpml) | 场更新主循环内**内嵌**记忆变量 ψ 的递推与修正项(Ex/Ey/Ez 与 Hx/Hy 共 24 个 ψ);Hz 的 z 方向采用基于 `c_h_zz` 的递归卷积(非 ψ)    |
-
-### 4.2 参数与含义
-
-| 参数                | 默认值                        | 含义                                       |
-| ------------------- | ----------------------------- | ------------------------------------------ |
-| `PML_X,PML_Y,PML_Z` | input.dat 第 5 行(建议 ≥5 层) | 三个方向的 PML 层数                        |
-| `ma`                | 3                             | σ 沿厚度方向的多项式阶数(由内向外幂律增长) |
-| `mb`                | 1                             | α 沿厚度方向的多项式阶数                   |
-| `sig_max`           | 1.0e2                         | PML 外侧最大电导率(决定吸收强度)           |
-| `alpha_max`         | **1.0e-1**                    | 复频移因子最大值                           |
-| `kappa_max`         | 1.0                           | 坐标拉伸系数最大值(1 表示不拉伸)           |
-
-σ、α 沿厚度的空间分布(以 x 方向下层为例,其余边界对称):
-
-```
-σ(i) = sig_max · ((L−i)/(L−1))^ma
-α(i) = alpha_max · ((i−1)/(L−1))^mb
-```
-
-E 场采样在整层、H 场采样在半层(交错),因此 H 方向的 σ/α/κ 按半层偏移
-(`i−0.5`)构造,与 E 方向错开。
-
-### 4.3 使用说明
-
-1. `input.dat` 第 4 行 `Logic_PML=1`,第 5 行给出三个方向的 PML 层数
-   (如 `15,15,15`);PML 层内网格尺寸应与核心区一致(保持均匀)。
-2. **α_max 取 0.1 是关键调参**:α(复频移因子)负责吸收低频扩散场。
-   若取值过小(如 0.01),晚时段的低频反射场不能及时衰减,会在边界往返
-   叠加,导致关断后约 10⁻⁵ s 量级出现指数发散(结果为 NaN)——这是
-   CPML 版最常见的不稳定来源,务必保持 α_max=1.0e-1。
-3. 切换回原始边界:第 4 行改为 `0` 即可,行为与旧版本完全一致,无需重新编译。
-
-### 4.4 与 Dirichlet 边界的对比
-
-|          | CPML(`Logic_PML=1`)               | Dirichlet(`Logic_PML=0`) |
-| -------- | --------------------------------- | ------------------------ |
-| 边界处理 | 吸收层,模拟无界空间               | 边界处场直接为零         |
-| 晚时精度 | 吸收反射,衰减曲线平直             | 边界反射可能污染晚时响应 |
-| 计算量   | 每步多 24 个记忆变量递推(约 +30%) | 无额外开销               |
-| 网格要求 | PML 层内需均匀网格                | 无特殊要求               |
-| 稳定性   | 调参正确时稳定                    | 稳定                     |
-
-> **一致性验证**:均匀半空间算例(81×81×80 网格、PML 15 层、1 ms 平台、
-> 关断后对比)中,CPML 与 Dirichlet 两种边界在关断后早期(场尚未到达边界
-> 时)的响应曲线一致,差异 <0.2%(源自两版循环次序不同导致的浮点舍入累积,
-> 非物理差异),证明 CPML 实现与主迭代等价、正确。
-
---
-
-## 5. 输出文件说明
-
-| 文件                            | 内容                                                                                                              |
-| ------------------------------- | ----------------------------------------------------------------------------------------------------------------- |
-| `dBzdt_1.txt`, `dBzdt_2.txt`, … | 每个接收点一个文件。文件头 2 行为说明(测点编号、测点坐标),其后每行 3 列:迭代步数、关断后时间(s)、该时刻磁场响应值 |
-| `CTIME_TIXING_UPCOS.DAT`        | 计算时间序列,每行 3 列:累计时间、时间步长、源电流幅值(波形)                                                       |
-| `conductivity.vtk`              | 模型电导率分布(规则网格 VTK 格式),可用 ParaView/Tecplot 等打开,检查模型是否正确构建                               |
-| `logfile.log`                   | 运行日志:参数回显、格式选择提示、运行错误等                                                                       |
-| `fort.5141`                     | 共形网格计算过程的调试输出                                                                                        |
-| `TEM_decay_curve.png`           | 衰减曲线图(由第 [6 节](#6-衰减曲线快速成图tem_decay_plotpy)的 `TEM_decay_plot.py` 生成)                           |
-
-运行结束时屏幕会打印各分段的迭代进度、总计算耗时;正常完成后 `logfile.log`
-末尾出现 `Computation finished!`。
-
---
-
-## 6. 衰减曲线快速成图(TEM_decay_plot.py)
-
-程序目录下的 `TEM_decay_plot.py` 用于将正演结果 `dBzdt_*.txt` 快速绘制为
-**衰减电压曲线图**(双对数坐标)。
-
-### 6.1 使用方式
-
-```bash
-# 需要 numpy 与 matplotlib
-pip install numpy matplotlib
-
-# 在计算输出文件(dBzdt_*.txt)所在目录运行
-python TEM_decay_plot.py
-```
-
-脚本自动搜索脚本同目录下所有 `dBzdt_*.txt` 文件,每个接收点画一条曲线,
-测点编号与坐标自动标注在图例中;默认输出高分辨率图片 **`TEM_decay_curve.png`**
-(dpi=600)并弹窗显示。
-
-### 6.2 主要可调参数(脚本顶部"User parameters"区)
-
-| 参数                | 默认值                | 说明                                                             |
-| ------------------- | --------------------- | ---------------------------------------------------------------- |
-| `file_pattern`      | `dBzdt_*.txt`         | 匹配的结果文件模式                                               |
-| `xmin, xmax`        | `1e-6, 1e-1`          | 横轴(时间, s)显示范围                                            |
-| `ymin, ymax`        | `None, None`          | 纵轴(响应)显示范围,`None` 表示自动                               |
-| `use_abs`           | `True`                | `True` 画 \|dBz/dt\|(正响应),`False` 画带符号值;纵轴标签随之切换 |
-| `savefig`           | `True`                | 是否保存图片                                                     |
-| `save_name`         | `TEM_decay_curve.png` | 保存文件名                                                       |
-| `dpi`               | `600`                 | 图片分辨率                                                       |
-| `label_fontsize` 等 | 18/15/15              | 标签、刻度、图例字号                                             |
-| `linewidth`         | `2.5`                 | 曲线线宽                                                         |
-
-超过 8 条曲线时自动改用顺序蓝色渐变配色,不会循环重复颜色。
-
-> 曲线取值规则:仅绘制 `时间 > 0` 且响应值 > 0 的点(双对数坐标下负值/零值无法显示)。
-
---
-
-## 7. 快速上手(本目录自带算例)
-
-本目录自带算例:**起伏地形下的复杂三维异常体**(模型 101×101×100 网格、网格尺寸
-20 m、源边长 500 m、背景 0.01 S/m、异常体为低阻体 4.0 S/m、1 个测点(源中心正
-下方);异常体网格 2663 节点/5322 单元,地形网格 10039 节点/5426 单元)。
-
-运行步骤:
-
-1. 确认目录下存在:`input.dat`、`Complex_anomalous.dat`(及 `.stl`,描述同一异常体)、
-   `Complex_Terrain.dat`(及 `.stl`,描述同一地形)。两类文件 `.dat` 与 `.stl`
-   同时存在时,程序以 `.dat` 优先并给出提示。
-2. VS2019 打开 `tem3dfdtd-cmm.sln` → 选择 `Release | x64` → 生成。
-3. 将生成的 exe 复制到本目录(或把输入文件放入 exe 目录)后运行。
-4. 观察屏幕输出,正常流程为:
-   ```
-   Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.
-   Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.
-   The number of grids in the core area is odd
-   At least 320M memory is needed!
-   ...
-   Conformal mesh of terrain is complete!
-   Ray tracing computation of terrain is complete!
-   Conformal mesh of terrain is finished
-   Start conformal processing of the anomalous body
-   ...
-   Now computing fraction:  1
-   50 steps have just finished
-   ...
-   ```
-5. 计算完成后检查输出文件 `dBzdt_1.txt` 、`conductivity.vtk`。
-6. (可选)运行 `python TEM_decay_plot.py` 生成衰减曲线图 `TEM_decay_curve.png`
-   (见第 [6 节](#6-衰减曲线快速成图tem_decay_plotpy))。
-
---
-
-## 8. 常见问题
-
-**Q1:运行提示 `libiomp5md.dll` 找不到**
-OpenMP 运行时库缺失。将 Intel oneAPI 安装目录 `bin/libiomp5md.dll` 复制到 exe 旁
-(或加入 PATH)。
-
-**Q2:提示 `Both Complex_Terrain.dat and Complex_Terrain.stl exist! ...`**
-两个格式文件都在。程序以 `.dat` 优先。若想用 STL,请将 `.dat` 文件移走或改名。
-
-**Q3:计算很慢 / 内存不足**
-减少 `NX,NY,NZ` 或增大 `GridSize`;控制 `NSTOP`;`MAX_OFF_TIME` 决定实际迭代步数,
-程序会以两者中的较小者为准。运行前会打印所需内存估算。
-
-**Q4:如何只算均匀半空间(无异常体、无地形)?**
-半空间模型应包含"空气 + 大地"两部分。将 `Complex_anomalous.*` 与
-`Complex_Terrain.*` 移走,并在 `input.dat` 中设 `TEMP_II = 2`:第 1 块设为上半部分
-(空气,电导率如 `1e-5`),第 2 块设为下半部分(大地,电导率如 `0.01`),即构成
-均匀半空间。注意:**`TEMP_II = 0` 时整个模型只填充背景电导率(全空间均匀介质,
-不含空气层)**。
-
-**Q5:接收点坐标怎么写?**
-坐标是**相对回线源中心**的局部坐标(单位 m),正负方向与坐标轴一致。
-
-**Q6:CPML 吸收边界与原始 Dirichlet 边界怎么选?**
-`input.dat` 第 4 行开关 `Logic_PML`:`1` 启用 CPML 吸收边界(第 5 行
-`10,10,10` 为三个方向的 PML 层数,可自行调整),能有效吸收边界反射,晚时
-(大偏移/晚时间)衰减曲线更平直;`0` 使用原始非均匀网格 Dirichlet(零场)
-边界。切换开关无需重新编译。
-
---
-
-## 9. 参考文献
-
-[1] 孙怀凤, 李貅, 李术才, 等. 考虑关断时间的回线源激发TEM三维时域有限差分正演[J]. 地球物理学报, 2013, 56(3): 1049-1064.
-
-[2] 柳尚斌, 李雪峰, 蓝日彦, 等. 瞬变电磁低频近似Maxwell方程的CPML吸收边界及施加方法[J]. 地球物理学报, 2022, 65(4): 1472-1481.
-
-[3] Li X, Zhao Q, Hu S, et al. Introducing complex geometries to Yee cells in FDTD for transient electromagnetic forward modeling[J]. Geophysics, 2025, 91(2): F1-F12.
@@ -1,16 +1,18 @@
 this folder is used to store the examples!

-专门为开源地球物理学校准备的算例：
-ex1 均匀半空间模型
-ex2 均匀半空间中包含一个低阻球体模型
-ex3 起伏地形下包含复杂异常体的模型
+# 专门为开源地球物理学校准备的算例：
+
+- ex1 均匀半空间模型
+- ex2 均匀半空间中包含一个低阻球体模型
+- ex3 起伏地形下包含复杂异常体的模型
+
+## 在每个算例文件夹下的：

-在每个算例文件夹下的：
 | 文件                                        | 内容                                                                                                                                                                 |
-| ------------------------------- | ----------------------------------------------------------------------------------------------------------------- |
-| input???.dat | 代表输入文件，其中Dirichlet代表采用Dirichlet边界条件，PML代表采用PML吸收边界条件 |
+| ------------------------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| input???.dat                                | 代表输入文件，其中Dirichlet代表采用Dirichlet边界条件，PML代表采用PML吸收边界条件，运行时将输入文件放入**tem3dfdtd-open\tem3dfdtd**目录中，并将文件名修改为 input.dat |
 | Complex_anomalous.dat Complex_anomalous.stl | 代表复杂异常体文件                                                                                                                                                   |
 | Complex_Terrain.txt Complex_Terrain.stl     | 代表复杂地形体文件                                                                                                                                                   |
 | 文件夹                                      | 代表计算结果                                                                                                                                                         |

-请参考使用
+> 所有算例请参考使用
@@ -1,21 +1,21 @@
 1			! 1-Ground TEM；2-SATEM
-500			! 回线源尺寸
-61,61,60	! x、y、z网格数量
-1			! 边界条件开关：1=CPML吸收边界，0=原始非均匀网格Dirichlet边界
+210			! 回线源尺寸
+121,121,120	! x、y、z网格数量
+0			! 边界条件开关：1=CPML吸收边界，0=原始非均匀网格Dirichlet边界
 10,10,10		! PML层数（x、y、z方向，仅开关=1时有效）
-25,25		! x方向均匀核心区域范围/数量
-25,25		! y方向均匀核心区域范围/数量
-25,25		! z方向均匀核心区域范围/数量
-20			! 核心均匀网格尺寸
+30,30		! x方向均匀核心区域范围/数量
+30,30		! y方向均匀核心区域范围/数量
+20,30		! z方向均匀核心区域范围/数量
+10			! 核心均匀网格尺寸
 0.01			! 背景介质电导率(S/m)
 2			! 介质数量
-1,61		
-1,61			
-1,30			
+1,121		
+1,121			
+1,60			
 1e-5	
-1,61
-1,61
-31,60
+1,121
+1,121
+61,120
 1e-2
 4000000		! 最大迭代次数
 90.101		! 最大计算时间
@@ -26,8 +26,10 @@
 1   			! AMP
 2.0 			! The conductivity of abnormal body (When the file Complex_anomalous.dat exists)
 TIXING_UPCOS    ! source type
-2                      ! 接收测点数量   
-1
-0,0,0
+3                      ! 接收测点数量   
+1			! 测点1编号
+-60,0,0			! 测点1坐标
 2
-0,0,-20
+-30,0,0
+3
+0,0,0
@@ -1,64 +0,0 @@
-!the standard input file for TEM calculation.Version 2.0 start to creat @2016-10-30 by Huaifeng Sun
-!finishe to creat @
-
-
-!this is the configuration type for the modeling
-!airborne, semi-airborne, ground, surface-borehole, tunnel, marine are the options
-#configuration_type:
-SEMI
-
-!source length
-#source-parameters:
-3
-
-!cell numbers in x,y,z directions
-!minumun grid size in uniform parts
-#cell-parameters:
-161,161,160
-0.5
-
-!background resistivity
-!number if abnormal bodies
-!the numbered abnormal body cell range in x diretion,y diretion,z diretion,resistivity
-#resistivity-parameters:
-0.01
-6
-1,161,1,161,1,80,1e-4
-1,161,1,161,80,81,0.002
-1,161,1,161,81,83,0.0033
-76,86,76,86,86,92,0.02
-76,86,76,86,92,98,0.033
-76,86,76,86,98,106,0.05
-
-#n-stop:
-3500000
-
-!The maximum computation time, unit of which is ms
-!the raise time and its step
-!the wave length time
-!the ramp time and its step
-!timestep
-!current in amper
-!waveform-type
-#waveform-parameters:
-30.002
-1e-6,1e-9
-10000e-6
-1e-6,1e-9
-1e-7
-1
-TIXING_UPCOS
-
-!number of flight hight
-!flight hight 
-#flight-parameters:
-2
-0.5,1
-
-!HE stands for horizontal value,HZ stands for vertical value
-!number of survey lines
-!start and stop point cell number
-#receiver-parameters:
-HE
-1
-66,96
@@ -1,43 +0,0 @@
-SEMI
-3
-161,161,160
-0.5
-0.01
-6
-1,161
-1,161
-1,80
-1e-4
-1,161
-1,161
-80,81
-0.002
-1,161
-1,161
-81,83
-0.0033
-76,86
-76,86
-86,92
-0.02
-76,86
-76,86
-92,98
-0.033
-76,86
-76,86
-98,106
-0.05
-3500000
-30.002
-1e-6,1e-9
-10000e-6
-1e-6,1e-9
-1e-7
-1
-2
-0.5,1
-TIXING_UPCOS
-HE
-1
-66,96
@@ -1,99 +0,0 @@
-<?xml version="1.0" encoding="utf-8"?>
-<!-- Version 1.0-->
-<model>
-  
-  <!--!this is the configuration type for the modeling-->
-  <!--!airborne, semi-airborne, ground, surface-borehole,
-    tunnel, marine are the options-->
-  <configuration>SEMI</configuration>
-
-  <source>
-    <source_length>3</source_length>
-    <source_type>loop</source_type>
-    <current>1.0</current>
-    <Tx_waveform>TIXING_UPCOS</Tx_waveform>
-  </source>
-  <cell>
-    <x>161</x>
-    <y>161</y>
-    <z>160</z>
-    <mini_size>0.5</mini_size>
-  </cell>
-
-  <!--This is the time parameters in the iteration. -->
-  <!--the total_time corresponds to the computation time from 0 to the end-->
-  
-  <time>
-    <total_time>30.002</total_time >
-    <raise_time>1e-6,1e-9</raise_time >
-    <ramp_time>1e-6,1e-9</ramp_time>
-    <on_time>1.0e-2</on_time>
-    <off_time_step_max>1.0e-7</off_time_step_max>
-    <max_stop_iteration>3500000</max_stop_iteration>
-  </time>
-
-  <Rx_parameters>
-    <num_record_height>2</num_record_height>
-    <heights>0.5,1</heights>
-    <rec_flag>HE</rec_flag>
-    <Rx_lines_num>1</Rx_lines_num>
-    <Rx_range>66,96</Rx_range>
-  </Rx_parameters> 
-  
-  <resistivity_parameters>
-    <background>0.01</background>
-    <abnormal>
-      <number>6</number>
-      <block1>
-        <x_range>1,161</x_range>
-        <y_range>1,161</y_range>
-        <z_range>1,80</z_range>
-        <conductivity>1e-4</conductivity>
-      </block1>
-
-      <block2>
-        <x_range>1,161</x_range>
-        <y_range>1,161</y_range>
-        <z_range>80,81</z_range>
-        <conductivity>0.002</conductivity>
-      </block2>
-
-      <block3>
-        <x_range>1,161</x_range>
-        <y_range>1,161</y_range>
-        <z_range>81,83</z_range>
-        <conductivity>0.0033</conductivity>
-      </block3>
-
-      <block4>
-        <x_range>76,86</x_range>
-        <y_range>76,86</y_range>
-        <z_range>86,92</z_range>
-        <conductivity>0.02</conductivity>
-      </block4>
-
-      <block5>
-        <x_range>76,86</x_range>
-        <y_range>76,86</y_range>
-        <z_range>82,98</z_range>
-        <conductivity>0.033</conductivity>
-      </block5>
-
-      <block6>
-        <x_range>76,86</x_range>
-        <y_range>76,86</y_range>
-        <z_range>98,106</z_range>
-        <conductivity>0.05</conductivity>
-      </block6>
-      
-    </abnormal>
-     
-    
-  </resistivity_parameters>
-  
-
-
-
-
-
-</model>
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
 SUBROUTINE anomalous_conformal
 !This procedure is used to compute anomalous body conformal
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 !function description
    !this suboutine is used to select and call some suboutine to close some 
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine GetSourcePosition
  ! This subroutine initialize the value of array is_ex_in_source which is used when judging whether the grid contains the source
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Shangbin Liu(lsbin87@126.com)
    
 SUBROUTINE Get_Receiver_Gridlabel
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 ! --------------------------------Subroutine part---------------------------------------------!
 subroutine Iteration_cpml
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine OpenRecFiles
  ! This subroutine opens all the files needed to record data of interests
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine SubCloseRecFiles(Flag)
  use constantparameters
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine SubOpenRecFiles(Flag)
  use constantparameters
@@ -72,7 +72,7 @@ subroutine SubOpenRecFiles(Flag)
  end select
  ! ------------------------------------end of File name Distribution---------------------------------------------!
  ! --------------------------------------open file code-------------------------------------------!
-  ! if the compiler reports the error: 'Too Many Open FIles!', you can come to tdem.org website and find the solutions.
+  ! if the compiler reports the error: 'Too Many Open FIles!', you can come to https://git.em3d.cn/ website and find the solutions.
  ! -----------------------------------------------------------------------------------------------------!
  do ii=1,NumRecHeights+1,1
     do jj=1,NumRecLines,1
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine SubWriteRecFiles(Flag,num)
  ! This subroutine writes all the data of intersted recording points from different recording plane which are given in the input.dat file.
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
 SUBROUTINE terrain_conformal

@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 !function description
    !this subroutine is used to allocate dynamic memory to the selected array.
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 !function description
    !this suboutine is used to write out the readed calculation parameters for 
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 !function description
    !this suboutine is used to close #5300 file.
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE FREE_MEMORY
  USE CONSTANTPARAMETERS
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 subroutine Get_eps_r
  use constantparameters
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine Get_mstop
  ! if the value of plus is too small, it will cause array bounds exceeded because the array bonds of Mstop and Mstart is set to 10000.
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 subroutine Get_pml_parameters
  use constantparameters
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 !SUBROUTINE GET_SYS_TIMEDATA(OUTPUT)
 !  ! The original subroutien written by Huaifeng Sun can not work in PGI compiler, so I change it into what it looks like here.
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 SUBROUTINE GET_COORDINATES
 	!> @brief This subroutine calculates the coordinates of each grid node in the 3D domain,
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE GET_NON_UNIFORMGRID
  USE CONSTANTPARAMETERS
@@ -30,17 +30,17 @@ SUBROUTINE GET_NON_UNIFORMGRID
  ! ----------------------------end of mesh------------------------------------!
  ! -------------------------------mesh x----------------------------------------!
  IF(Logi_Sourcelenth) THEN    !The number of grids occupied by the source is odd
-	if(SourceLength/GridSize.gt.51)then
-     do ii=nxs-(SourceLength/GridSize-1)/2,nxs+(SourceLength/GridSize-1)/2,1
+	if(SourceGridNum.gt.51)then
+     do ii=nxs-(SourceGridNum-1)/2,nxs+(SourceGridNum-1)/2,1
        Cdelx(ii)=GridSize
     end do
-     do ii=nxs-(SourceLength/GridSize-1)/2-1,1,-1
+     do ii=nxs-(SourceGridNum-1)/2-1,1,-1
        Cdelx(ii)=Cdelx(ii+1)*scale_par
        if(Cdelx(ii).gt.GridSize_MAX)then
           Cdelx(ii)=GridSize_MAX
        end if
     end do
-     do ii=nxs+(SourceLength/GridSize-1)/2+1,nx,1
+     do ii=nxs+(SourceGridNum-1)/2+1,nx,1
        Cdelx(ii)=Cdelx(ii-1)*scale_par
        if(Cdelx(ii).gt.GridSize_MAX)then
           Cdelx(ii)=GridSize_MAX
@@ -64,17 +64,17 @@ SUBROUTINE GET_NON_UNIFORMGRID
     end do
 	ENDIF
  ELSE    !The number of grids occupied by the source is even
-	if(SourceLength/GridSize.gt.51)then
-		do ii=nxs-(SourceLength/GridSize-1)/2,nxs+1+(SourceLength/GridSize-1)/2,1
+	if(SourceGridNum.gt.51)then
+		do ii=nxs-(SourceGridNum-1)/2,nxs+1+(SourceGridNum-1)/2,1
 			Cdelx(ii)=GridSize
 		end do
-		do ii=nxs-(SourceLength/GridSize-1)/2-1,1,-1
+		do ii=nxs-(SourceGridNum-1)/2-1,1,-1
 			Cdelx(ii)=Cdelx(ii+1)*scale_par
 			if(Cdelx(ii).gt.GridSize_MAX)then
 				Cdelx(ii)=GridSize_MAX
 			end if
 		end do
-		do ii=nxs+(SourceLength/GridSize-1)/2+2,nx,1
+		do ii=nxs+(SourceGridNum-1)/2+2,nx,1
 			Cdelx(ii)=Cdelx(ii-1)*scale_par
 			if(Cdelx(ii).gt.GridSize_MAX)then
 				Cdelx(ii)=GridSize_MAX
@@ -101,17 +101,17 @@ SUBROUTINE GET_NON_UNIFORMGRID
  ! -----------------------------end of mesh------------------------------------!
  ! --------------------------------mesh y----------------------------------------!
  IF(Logi_Sourcelenth) THEN    !The number of grids occupied by the source is odd
-	IF(SourceLength/GridSize.gt.51)then
-		do ii=nys-(SourceLength/GridSize-1)/2,nys+(SourceLength/GridSize-1)/2,1
+	IF(SourceGridNum.gt.51)then
+		do ii=nys-(SourceGridNum-1)/2,nys+(SourceGridNum-1)/2,1
 			Cdely(ii)=GridSize
 		enddo
-		do ii=nys-(SourceLength/GridSize-1)/2-1,1,-1
+		do ii=nys-(SourceGridNum-1)/2-1,1,-1
 			Cdely(ii)=Cdely(ii+1)*scale_par
 			if(Cdely(ii).gt.GridSize_MAX)then
 				Cdely(ii)=GridSize_MAX
 			end if
 		end do
-		do ii=nys+(SourceLength/GridSize-1)/2+1,ny,1
+		do ii=nys+(SourceGridNum-1)/2+1,ny,1
 			Cdely(ii)=Cdely(ii-1)*scale_par
 			if(Cdely(ii).gt.GridSize_MAX)then
 				Cdely(ii)=GridSize_MAX
@@ -135,17 +135,17 @@ SUBROUTINE GET_NON_UNIFORMGRID
 		end do
 	end if
  ELSE	!The number of grids occupied by the source is even
-	IF(SourceLength/GridSize.gt.51)then
-		do ii=nys-(SourceLength/GridSize-1)/2,nys+1+(SourceLength/GridSize-1)/2,1
+	IF(SourceGridNum.gt.51)then
+		do ii=nys-(SourceGridNum-1)/2,nys+1+(SourceGridNum-1)/2,1
 			Cdely(ii)=GridSize
 		enddo
-		do ii=nys-(SourceLength/GridSize-1)/2-1,1,-1
+		do ii=nys-(SourceGridNum-1)/2-1,1,-1
 			Cdely(ii)=Cdely(ii+1)*scale_par
 			if(Cdely(ii).gt.GridSize_MAX)then
 				Cdely(ii)=GridSize_MAX
 			end if
 		end do
-		do ii=nys+(SourceLength/GridSize-1)/2+2,ny,1
+		do ii=nys+(SourceGridNum-1)/2+2,ny,1
 			Cdely(ii)=Cdely(ii-1)*scale_par
 			if(Cdely(ii).gt.GridSize_MAX)then
 				Cdely(ii)=GridSize_MAX
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 SUBROUTINE GETDATA
  USE CONSTANTPARAMETERS
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
 !this subroutine is writen by Huaifeng Sun from May 29, 2017

 subroutine getxmldata
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE MEMORY_USE_ESTIMATION 
  !This subroutine is written by Huaifeng Sun, and it has not been modified since the last ice age, so it can not reveal the real consumption of memory now.
@@ -11,17 +11,17 @@ SUBROUTINE MEMORY_USE_ESTIMATION
  CHARACTER(LEN=40) XSTRING
  INTEGER(KIND=1) CONTD
  TUSE=0.0
-  TUSE=TUSE+NX*NYB*NZB+NXB*NY*NZB+NXB*NYB*NZ      !电场E使用内存
-  TUSE=TUSE+NXB*NY*(NZ+1)+NX*NYB*(NZ+1)+NX*NY*NZB !磁场H使用内存
-  TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !边界条件使用内存
-  TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !边界条件使用内存
-  TUSE=TUSE+NXB*NYB*NZB   !模型使用内存
+  TUSE=TUSE+NX*NYB*NZB+NXB*NY*NZB+NXB*NYB*NZ      !�糡Eʹ���ڴ�
+  TUSE=TUSE+NXB*NY*(NZ+1)+NX*NYB*(NZ+1)+NX*NY*NZB !�ų�Hʹ���ڴ�
+  TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !�߽�����ʹ���ڴ�
+  TUSE=TUSE+4*NY*NZB+4*NYB*NZ+NX*4*NZB+NXB*4*NZ+NX*NYB*4+NXB*NY*4 !�߽�����ʹ���ڴ�
+  TUSE=TUSE+NXB*NYB*NZB   !ģ��ʹ���ڴ�
  TUSE=TUSE+NSTOP*2
  TUSE=TUSE/1024
  TUSE=TUSE/1024
  TUSE=TUSE*16
  WRITE (XSTRING,'(I40)') TUSE
-  XSTRING = 'At least '//TRIM(ADJUSTL(XSTRING))//'M memory is needed!' !拼接为要求的FORMAT格式
+  XSTRING = 'At least '//TRIM(ADJUSTL(XSTRING))//'M memory is needed!' !ƴ��ΪҪ���FORMAT��ʽ
  XSTRING = TRIM(ADJUSTL(XSTRING)) 
  WRITE(*,*)XSTRING
  RETURN
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
    
 SUBROUTINE RES_CONFIGURE
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    

 SUBROUTINE SIN_SOURCE
@@ -25,6 +25,6 @@ SUBROUTINE SIN_SOURCE
     WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I)
  ENDDO
  CLOSE(9)
-  WRITE(10005,*)'半正弦发射波形时间序列已经写入文件CTIME_SIN_SOURCE.DAT'
+  WRITE(10005,*)'�����ҷ��䲨��ʱ�������Ѿ�д���ļ�CTIME_SIN_SOURCE.DAT'
  RETURN
 ENDSUBROUTINE SIN_SOURCE
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE TIME_SERIOUS
  USE CONSTANTPARAMETERS
@@ -8,7 +8,7 @@ SUBROUTINE TIME_SERIOUS
  USE RES_MODEL_PARAMETER
  USE TIME_PARAMETER
  USE OMP_LIB
-  !本子程序将计算初始时间和时间序列
+  !���ӳ��򽫼����ʼʱ���ʱ������
  IMPLICIT NONE
  INTEGER NSTOP_TEMP,i,j,k
  TIME_MAX=100*GridSize*SQRT(EPS0*MU0/3.0) !Time_max can be set to larger value if the value of GridSize if less than 1m, otherwise you will spend a much longer time in calculation.
@@ -42,27 +42,27 @@ SUBROUTINE TIME_SERIOUS
  ENDDO                         !This subroutine computes the value of Nstop which satisfies the requirement of Max_off_time
  IF(NSTOP_TEMP .LT. NSTOP)THEN
     NSTOP=NSTOP_TEMP           !Change the value of Nstop to a smaller value according to the above computation
-     WRITE(10005,*)'NSTOP改变为',NSTOP        
+     WRITE(10005,*)'NSTOP�ı�Ϊ',NSTOP        
     OPEN(9,FILE='CTIME_TIXING_UPCOS.DAT',STATUS='UNKNOWN')
     DO I=1,NSTOP
        WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I)
     ENDDO
     CLOSE(9)
-     WRITE(10005,*)'采用升余弦函数和降余弦函数处理的梯形波发射波形时间序列已经写入文件CTIME_TIXING_UPCOS.DAT'
+     WRITE(10005,*)'���������Һ����ͽ����Һ������������β����䲨��ʱ�������Ѿ�д���ļ�CTIME_TIXING_UPCOS.DAT'
  ELSEIF(NSTOP_TEMP .EQ. NSTOP)THEN
     NSTOP=NSTOP_TEMP
-     WRITE(10005,*)'NSTOP没有改变，可能无法满足最大最大延时设置，请重新设置NSTOP.'
+     WRITE(10005,*)'NSTOPû�иı䣬�����޷�������������ʱ���ã�����������NSTOP.'
     OPEN(9,FILE='CTIME_TIXING_UPCOS.DAT',STATUS='UNKNOWN')
     DO I=1,NSTOP
        WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I)
     ENDDO
     CLOSE(9)
-     WRITE(10005,*)'采用升余弦函数和降余弦函数处理的梯形波发射波形时间序列已经写入文件CTIME_TIXING_UPCOS.DAT'
+     WRITE(10005,*)'���������Һ����ͽ����Һ������������β����䲨��ʱ�������Ѿ�д���ļ�CTIME_TIXING_UPCOS.DAT'
  ELSE
     WRITE(10005,*)'The number of iteration steps exceeds the range given in the input.dat, please change it. now the Nstop value is determined by Max_off_time.'
     print*,'The number of iteration steps exceeds the range given in the input.dat, please change it. now the Nstop value is determined by Max_off_time.'
     print*,'I give you a pause here, you should decide to continue or to quit'
-     pause
+     READ(*,*)
     Nstop=Nstop_temp
  ENDIF
  RETURN
@@ -1,8 +1,8 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
-!梯形波电流激发，采用升余弦开关函数和降余弦开关函数
+!���β��������������������ҿ��غ����ͽ����ҿ��غ���
 SUBROUTINE TIXING_SOURCE_UPCOS
  USE CONSTANTPARAMETERS
  USE ELECTROMAGNETIC_VARIABLES
@@ -20,16 +20,16 @@ SUBROUTINE TIXING_SOURCE_UPCOS
  SOURCE(1)=AMP*0.5*(1-COS(PI*CTIME(1)/RAISETIME))  !AMP*CTIME(1)/RAISETIME
  DO I=2,NSTOP
     CTIME(I)=CTIME(I-1)+DELT(I-1)
-     IF(CTIME(I) .LT. RAISETIME)THEN !小于raistime
+     IF(CTIME(I) .LT. RAISETIME)THEN !С��raistime
        DELT(I)=RAISESTEP
        SOURCE(I)=AMP*0.5*(1-COS(PI*CTIME(I)/RAISETIME))  !AMP*CTIME(I)/RAISETIME
-     ELSEIF(CTIME(I) .GE. RAISETIME .AND. CTIME(I) .LT. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM)THEN  !RAISETIME小于等于CTIME(I)<RAISETIME+WAVE-TIME_RAMP2WAVE_SUM
+     ELSEIF(CTIME(I) .GE. RAISETIME .AND. CTIME(I) .LT. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM)THEN  !RAISETIMEС�ڵ���CTIME(I)<RAISETIME+WAVE-TIME_RAMP2WAVE_SUM
        DELT(I)=DELT(I-1)*1.0005
        IF(DELT(I) .GE. WAVESTEP)THEN
           DELT(I)=WAVESTEP
        ENDIF
        SOURCE(I)=AMP  !1.0D0
-     ELSEIF(CTIME(I) .GE. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM .AND. CTIME(I) .LT. RAISETIME+WAVE)THEN !RAISETIME+WAVE-TIME_RAMP2WAVE_SUM小于等于CTIME(I)<RAISETIME+WAVE
+     ELSEIF(CTIME(I) .GE. RAISETIME+WAVE-TIME_RAMP2WAVE_SUM .AND. CTIME(I) .LT. RAISETIME+WAVE)THEN !RAISETIME+WAVE-TIME_RAMP2WAVE_SUMС�ڵ���CTIME(I)<RAISETIME+WAVE
        DELT(I)=DELT(I-1)*0.9995
        IF(DELT(I) .LE. RAMPSTEP)THEN
           DELT(I)=RAMPSTEP
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE TIXING_SOURCE
  USE CONSTANTPARAMETERS
@@ -47,6 +47,6 @@ SUBROUTINE TIXING_SOURCE
     WRITE(9,'(3E24.16E3)')CTIME(I),DELT(I),SOURCE(I)
  ENDDO
  CLOSE(9)
-  WRITE(10005,*)'梯形波发射波形时间序列已经写入文件CTIME_TIXING.DAT'
+  WRITE(10005,*)'���β����䲨��ʱ�������Ѿ�д���ļ�CTIME_TIXING.DAT'
  RETURN
 ENDSUBROUTINE TIXING_SOURCE
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE TRIANGLE_SOURCE
  USE CONSTANTPARAMETERS
@@ -26,6 +26,6 @@ SUBROUTINE TRIANGLE_SOURCE
     WRITE(9,'(3E24.16)')CTIME(I),DELT(I),SOURCE(I)
  ENDDO
  CLOSE(9)
-  WRITE(10005,*)'三角发射波形时间序列已经写入文件CTIME_TRIANGLE_SOURCE.DAT'
+  WRITE(10005,*)'���Ƿ��䲨��ʱ�������Ѿ�д���ļ�CTIME_TRIANGLE_SOURCE.DAT'
  RETURN
 ENDSUBROUTINE TRIANGLE_SOURCE
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 subroutine WriteRecFiles(num)
  use constantparameters
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 SUBROUTINE ZERO
  USE CONSTANTPARAMETERS
@@ -1,13 +1,13 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 ! This is a finite difference time domain (FDTD) code for the simulation of transient electromagnetic (TEM);
 ! This code is designed to  be used in semi_airborne TEM with a loop source;
-! This code is written by Huaifeng Sun (sunhuaifeng@gmail.com) and Xushan Lu (luxushan@gmail.com), The conformal mesh part is written by Xinyu Li (202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com);
+! This code is written by Huaifeng Sun (sunhuaifeng@email.sdu.edu.cn) and Xushan Lu (luxushan@gmail.com), The conformal mesh part is written by Xinyu Li (202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com);
 ! OpenACC API is used in this code for the acceleration with GPU device; therefore, you are recommended to compile this code with --
 ! --PGI Accelerator Fortran Workstation compiler. A Nvidia GPU card with CUDA capability is required if you want to run this code in parallel mode.
-! Nobody is allowed to copy or distribute this code to people outside of TDEM.org group without the permission from Prof. Xiu Li (lixiu@chd.edu.cn)--
+! Nobody is allowed to copy or distribute this code to people outside of https://git.em3d.cn/ group without the permission from Prof. Xiu Li (lixiu@chd.edu.cn)--
 ! --or you will be 
 ! Contact the author for more detailed information.

@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Qi Zhao(zhaoqi_326326@163.com)
 MODULE Precision
    USE ISO_FORTRAN_ENV, ONLY : i1k => INT8, i2k => INT16, i4k => INT32, i8k => INT64, &
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Qi Zhao(zhaoqi_326326@163.com)
 MODULE vtk_fix_header
    USE Precision, ONLY : i4k
@@ -1,4 +1,4 @@
-!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@gmail.com)
+!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
 !written by Qi Zhao(zhaoqi_326326@163.com)

 ! --------------------------------Subroutine part---------------------------------------------!
@@ -204,7 +204,7 @@ Module VTK_Fortran
                    do i=1,NXB
                        
                        temp = temp + 1
-                        write(unit, fmt='(3f)', iostat=error) Coord_x(i), Coord_y(j), Coord_z(k)
+                        write(unit, fmt='(3F25.16)', iostat=error) Coord_x(i), Coord_y(j), Coord_z(k)
                        
                    end do
                end do
@@ -250,7 +250,7 @@ Module VTK_Fortran
                do j=1, NY
                    do i=1,NX
                        temp = temp + 1
-                        write(unit, fmt='(9i)', iostat=error)nums,&
+                        write(unit, fmt='(9I12)', iostat=error)nums,&
                            ((i  ) + (j-1)*NXB + (k-1)*NXB*NYB) -1,  &
                            ((i+1) + (j-1)*NXB + (k-1)*NXB*NYB) -1,  &
                            ((i  ) + (j  )*NXB + (k-1)*NXB*NYB) -1,  &
@@ -266,7 +266,7 @@ Module VTK_Fortran
            write(unit=unit,fmt='(A,I0,A,I0)', iostat=error) "CELL_TYPES ", cells
            
            do i=1, temp
-                write(unit, fmt='(i)', iostat=error)num_cell_types
+                write(unit, fmt='(I12)', iostat=error)num_cell_types
            end do
            write(unit=unit,fmt='(2A,I0)', iostat=error) Value_site_cell," ", cells

@@ -311,7 +311,7 @@ Module VTK_Fortran
                do j=1, NY
                    do i=1,NX
                        temp = temp + 1
-                        write(unit, fmt='(f)', iostat=error)values(i,j,k)
+                        write(unit, fmt='(F25.16)', iostat=error)values(i,j,k)
                    end do
                end do
            end do
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
    MODULE CONSTANTPARAMETERS
  INTEGER NX,NY,NZ
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/

 MODULE ELECTROMAGNETIC_VARIABLES
  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE:: EX,EY,EZ    !The x,y and z component of electric field in 3 dimensions
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
 !The PML part follows the CPML scheme of Roden & Gedney (2000), integrated by the
 !author of the PML version (tem3dfdtd_第二版), with tuned parameters:
 !sig_x_max=1.0D2, alpha_x_max=1.0D-1, kappa_x_max=1.0 (RESTORED to the EXACT
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 MODULE RES_MODEL_PARAMETER
  REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE:: CCSIG,CCSIGX,CCSIGY,CCSIGZ    !The conductivity in each grid
@@ -1,6 +1,6 @@
-!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
-!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
-!Code distribution @ tdem.org or sunhuaifeng.com
+!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
+!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
+!Code distribution @ https://git.em3d.cn/
    
 MODULE TIME_PARAMETER
  REAL(KIND=8), DIMENSION(:), ALLOCATABLE:: CTIME,DELT
作者	SHA1	备注	提交日期
lixinyu	4c4b501612	修改了两处不严谨的地方，gfortran环境下会warning，1：get_non_uniformgrid.f90中，IF判断中IF(SourceLength/GridSize.gt.51)改SourceLength/GridSize为整型SourceGridNum；2：替换旧写法PAUSE语句为READ(,)	2026-08-18 18:48:15 +08:00
lixinyu	d855afd251	Merge branch 'master' of https://git.em3d.cn/openSource/tem3dfdtd-open	2026-08-18 16:21:44 +08:00
lixinyu	642849b079	修正一处bug，VTK_fortran.f90 格式化输出格式不规范，无法在gfortran环境中编译，现已改为标准格式	2026-08-18 16:19:47 +08:00
sunhf	d1d658598c	add English readme	2026-08-17 21:00:41 +08:00
sunhf	a1bfcdd26b	修正readme中的错误	2026-08-14 20:54:09 +08:00
sunhf	89ff871b0b	修正代码中的注释错误	2026-08-14 20:51:58 +08:00
lixinyu	5026369041	统一example测试案例中输入文件格式	2026-08-14 10:10:26 +08:00
sunhf	ceaf462d33	整合多余的readme	2026-08-14 08:05:08 +08:00