esa-sceva/satcom-chunk-collection huggingface dataset:
# A Study on Technology Competition of Graphene Biomedical Technology Based on Patent Analysis
## 4 Analysis  
### Patent Analysis  
#### 4.1.1 Annual Trend of Patent Activities  
The analysis of annual distribution of patent application and development trends of the top ten countries or regions leading the graphene biomedical race are presented in this section. Figure 2 illustrates the dynamic changes of patent application for graphene biomedical technology in China, the U.S., Korea, India, the U.K., Taiwan, Japan, Germany, Canada, France and the world. We find that the number of patents for graphene biomedical technology in the global had rapidly increased since 2012. In 2010, the Nobel Prize in physics was granted to [PERSON] and [PERSON] for their amazing breakthrough on graphene [59]. After that, many countries, with substantial efforts and investments from governments, have been focused towards the research and development (R&D) in graphene technology [60]. In addition, since the development trend of graphene has accelerated around the world, the number of graphene biomedical patents has rapidly grown. Regarding the number of patents, the top ten countries or regions hold 93.88% of global share in graphene biomedical patents. In addition, China (790 patents), the U.S. (378 patents), and Korea (271 patents) are the top three countries with more than 200 patents each. In total, the top three possess 83.08% of the global patent application in the graphene biomedical field. It is clear that China outstrips other countries or regions, and holds 45.61% of global share, suggesting its leading power in this field. One reason for this is the broad market and plentiful graphite resources in China, which offers great development opportunities for Chinese stakeholders [61]. The other reason is China's technology assessment system that put substantial effort into patents [61]. It is interesting to find that a country/region from Asia (e.g., China, Korea, Japan, and Taiwan) own the largest amount of graphene biomedical patents in total. Although the U.K, Germany and France play key roles in patenting activities in the European region, they are still far behind assignees from Asia and the U.S.  
In Figure 2, the U.S., Korea, and Japan started earlier in graphene biomedical technology, and filed their first patent in 2004. While in China, the patent application had appeared since 2005, followed by Taiwan in 2007, Germany and Canada in 2008. Before the year 2010, the U.S. (55%), Korea (12.5%), and Japan (12.5%) play the dominant roles in patenting activity, suggesting that these three countries obtain first-mover advantages. After 2011, the top three countries are China, the U.S., and Korea. China, with stable growth in the number of patents, becomes a technology catcher and takes the lion's share in the global market.
## 5 Experimental Results and Analysis
### Algorithm Timing Analysis  
In correlation operations, this paper utilizes the structure of FIR filters and optimizes them for multi-channel implementation. Within the FIR filter, the data are reused for computation, which avoids introducing delay but involves a trade-off between speed and area; thus, it is necessary to verify that the optimization does not affect the system's maximum clock frequency. This verification is conducted by analyzing the synthesis report of the hardware algorithm implementation. Table 4 presents the key parameter for the clock report, namely, the Worst Negative Slack (WNS). The results indicate that the _WNS_ values derived from the synthesis reports remain consistent at 3.313 ns before and after optimization. Using the formula \(F_{max}=1/(T-WNS)\), it is concluded that the maximum clock frequency \(F_{max}\) is not impacted by the FIR filter optimization, where \(T\) denotes the clock period.
# A Federated Generalized Linear Model for Privacy-Preserving Analysis
## 2 Materials and Methods
### Linear Regression
#### 2.3.2 Execution Time  
We were also interested in validating how the federated GLM's execution time escalates under different circumstances. For this purpose, we ran numerous simulations. All of them used the same mock server-client architecture and parameters described earlier.  
First, we explored the impact of increasing the number of _records_ while keeping the number of parties constant. Using the same corresponding data for each family as before, we simulated a _total_ of 30, 300, \(\ldots\), 3,000,000 records in a three-party scenario (i.e., each party had 10, 100, \(\ldots\), 1,000,000 records).  
Afterwards, we explored the impact of increasing the number of _parties_ while keeping the number of records constant. Using the same setup as before, we simulated scenarios with \(2,3,\ldots,10\) parties with a total of 10,000 records. Records were distributed randomly and practically evenly among parties (e.g., in a three-party scenario, two parties would have 3333 records, while the remaining one would have 3334).  
In both of these validations, in order to obtain a proper idea of the performance and to account for the variability due to inherent randomness in the data generation process, each case was simulated 100 times. We only measured the time (in seconds) from the beginning to the end of the federated GLM execution. In other words, we discarded the time it took to generate the data, since in a real-life situation, each party would already have their own data at hand.
# A Hybrid Deep Learning-Based Load Forecasting Model for Logical Range
## 3 Results
### Performance Evaluation Metrics  
To ensure a comprehensive and reliable evaluation of the GCSG model's performance, we adopt a multi-metric strategy comprising the coefficient of determination (\(R^{2}\), see Equation (12)), mean absolute error (MAE, see Equation (13)), and mean squared error (MSE, see Equation (14)). Specifically, \(R^{2}\) measures the proportion of variance explained by the model, while the MAE and MSE quantify the prediction errors from absolute and squared perspectives, respectively. This combination of metrics enables a more holistic assessment of the model's accuracy and robustness, avoiding overreliance on any single evaluation indicator.  
1. The coefficient of determination (\(R^{2}\)) measures the proportion of variance in the predicted Logical Range load data explained by the model. Higher values indicate better performance:  
\begin{table}
\begin{tabular}{l l} \hline
**Category** & **Specification** \\ \hline Operating system & Windows 10 (64-bit) \\ Programming language & Python 3.9 \\ Deep learning framework & PyTorch 2.0.0+ cu118 \\ Processor & 12 th Gen Intel(R) Core(TM) i9-12900K 3.20 GHz \\ GPU & NVIDIA GeForce RTX 3090 \\ Memory & 64 (GB) \\ \hline \end{tabular}
\end{table}
Table 1: Experimental setup.  
\[R^{2}=1-\frac{\sum\limits_{i=1}^{n}{(y_{i}-\hat{y}_{i})^{2}}}{\sum\limits_{i=1}^{n }{(y_{i}-\hat{y}_{i})^{2}}}. \tag{12}\]
2. The mean absolute error (MAE) is computed as the average absolute difference between predicted and actual load forecasting values: \[\text{MAE}=\frac{1}{n}\sum\limits_{i=1}^{n}{|y_{i}-\hat{y}_{i}|}.\] (13)
3. The mean squared error (MSE) quantifies average squared prediction errors, penalizing larger deviations: \[\text{MSE}=\frac{1}{n}\sum\limits_{i=1}^{n}{(y_{i}-\hat{y}_{i})^{2}}.\] (14)  
In the above formulas, \(\hat{y}_{i}\) denotes the predicted value, \(y_{i}\) denotes the true (actual) value, and \(\hat{y}\) is the mean of the true values.  
These metrics were used to assess the model's accuracy for time-series load forecasting in the Logical Range.
## 4 Scheme Design: MPMF
### Relay Selection
#### 4.2.1 Case 1
Here, for example, we have \(p_{1}=P_{(n_{1},d)}\), and \(p_{2}=P_{(n_{2},d)}\). If \(n_{x}\) and \(n_{y}\) have not contacted for a certain period, their DP values will decay and be updated by (4). \[P_{(x,y)}=P_{(x,y)old}\times\gamma^{k}, \tag{4}\]  
where \(\gamma\in[0,1]\) is an aging coefficient, and \(k\) denotes the number of expired time units after the last update of this DP value. Although the aging of \(P_{(x,y)}\) is ongoing after the last contact, \(P_{(x,y)}\) can still be increased by any other node \(n_{z}\) with \(P_{(x,z)}\) and \(P_{(z,y)}\) as:  
\[P_{(x,y)}=P_{(x,y)old}+(1-P_{(x,y)old})P_{(x,z)}P_{(z,y)}\times\beta. \tag{5}\]  
where \(\beta\) is a parameter with value \(\in[0,1]\). Let \(\overrightarrow{n_{i+1}}\) be a directional vector of \(n_{i+1}\), then a new \(\theta\) is calculated with respect to \(\overrightarrow{n_{i+1}}\) and \(\overrightarrow{n_{i+1}}d\). With \(\theta\), a secondary weighted form \(W^{(2)}\) is given as  
\[W^{(2)}_{i+1}=1-|\frac{\theta}{180}|. \tag{6}\]  
Then, combining (2) and (6), we have a two-hop weighed form \(W_{p}\) subject to a tuning parameter \(\alpha=(0,1)\) as  
Figure 2: Two-hop based routing between \(s\), \(n_{i}\), and \(n_{i+1}\). \[W_{p} =\alpha W_{i}^{(1)}+(1-\alpha)W_{i+1}^{(2)}\] \[=\alpha\times\frac{1}{\ln p_{i}\cdot\ln p_{i+1}+1}+(1-\alpha)\times 1 -|\frac{\theta}{180}|. \tag{7}\]  
If the resulting \(W_{p}\) is greater than a specific \(W_{0}\), \(n_{i}\) and \(n_{i+1}\) will be selected as candidate nodes to receive a message \(m_{i}\) from \(s\). First, \(n_{i}\) will receive \(m_{i}\) from \(s\) and then replicate it to \(n_{i+1}\). The above relay node selection can be written as follows:  
\[d(m_{i})= \begin{cases}1,&((\theta_{i}\wedge\theta_{i+1})\leq\frac{\pi}{2}) \wedge(\theta_{m}\leq\epsilon)\wedge(W_{p}>W_{0}),\\ 0,&\text{otherwise},\\ &\text{where }\epsilon\in[0,\frac{\pi}{2}],W_{0}\in[0,1],\text{ and }\theta_{m}=|\theta_{i}+\theta_{i+1}^{ \prime}|.\end{cases}  \tag{8}\]  
In (8), if both \(\theta_{i}\) and \(\theta_{i}\leq\frac{\pi}{2}\), \(\theta_{m}\leq\epsilon\), and \(W_{p}>W_{0}\), then both \(n_{i}\) and \(n_{i+1}\) will receive \(m_{i}\). In the case of \(W_{p}<W_{0}\), \(s\) will look for another one-hop neighbor node \(n_{i}^{\prime}\) and then repeat the above procedure to determine the next relay node. If no neighbor nodes can satisfy this selection policy, the original node keeps carrying the message during moving in OppNets.
## 4 Conclusions  
A study concerning a synthesis method based on a direct aerosolisation process of two solutions, each containing a reactant, has been proposed. The as-formed NPs, embedded in liquid droplets resulting from the coalescence of two different aerosolised phases, have been captured in a liquid phase by means of a Venturi nozzle. The present experimental technique relies on two closed circuits, namely a liquid phase devoted to NP abatement, and a gaseous phase, carrying two aerosolised liquid phases reacting with each other. The essentials of the present study can be summarised in the following points:  
* The use of two closed circuits, both gas- and liquid-tight, allows minimising NP dispersions in the surrounding environment. For this reason, this technique may be advantageously proposed in the synthesis of toxic or noxious NPs.
* As a consequence of the previous point, this synthesis process may be carried out in the presence of an inert gas, as a further possible extension of this method, without the need to replace gas leaks. For this reason, the present process may be economically attractive. Additionally, this aspect is of primary importance in the case of reactive NPs, like non-noble metal NPs, often requiring an inert atmosphere.
* The synthesis of NPs is carried out at room temperature, with positive effects in the synthesis of thermolabile NPs. No additional energy consumption and/or severe experimental conditions were needed, thus fulfilling inherent safety criteria for the process route and operator hazard exposure.
* This technique allowed obtaining inorganic NPs, avoiding contamination by surfactant or capping agents, usually adopted in wet chemical synthesis. For this reason, the present method lends itself to the production of NPs for theranostic uses, where the presence of additional compounds on the NPs' surfaces is often undesired.
* All other conditions being equal, especially in the absence of capping agents, the NPs synthesised here by direct aerosolisation proved to have a diameter dispersion narrower than those obtained by a wet chemical method based on liquid mixing.  
As a further development, the present technique will be tested in the synthesis of organic NPs, aiming to achieve inherent safer and sustainable application in nanomedicine.  
All authors contributed significantly and equally to all aspects pertaining to the realisation of the present paper. All authors have read and agreed to the published version of the manuscript.  
This research received no external funding.  
Not applicable.  
All data are reported in the paper.  
The authors wish to thank [PERSON] for the management of the use of the DLS instrument. The authors are grateful, for the FE-SEM characterizations, to the Laboratory of Electron Microscopy-LEM at the Department of Chemistry and Industrial Chemistry at the University of Genoa, Italy.  
The authors declare no conflicts of interest.
###### Abstract  
Underwater acoustic sensor networks play an important role in assisting humans to explore information under the sea. In this work, we consider the combination of sensor selection and data routing in three dimensional underwater wireless sensor networks based on Bayesian compressive sensing and particle swarm optimization. The algorithm we proposed is a two-tier PSO approach. In the first tier, a PSO-based clustering protocol is proposed to synthetically consider the energy consumption and uniformity of cluster head distribution. Then in the second tier, a PSO-based routing protocol is proposed to implement inner-cluster one-hop routing and outer-cluster multi-hop routing. The nodes selected to constitute \(i\)-th effective routing path decide which positions in the \(i\)-th row of the measurement matrix are nonzero. As a result, in this tier the protocol comprehensively considers energy efficiency, network balance and data recovery quality. The Bayesian Cramer-Rao Bound (BCRB) in such a case is analyzed and added in the fitness function to monitor the mean square error of the reconstructed signal. The experimental results validate that our algorithm maintains a longer life time and postpones the appearance of the first dead node while keeps the reconstruction error lower compared with the cutting-edge algorithms which are also based on distributed multi-hop compressive sensing approaches.  
Bayesian compressive sensing; particle swarm optimization; three dimensional underwater wireless sensor network; Bayesian Cramer-Rao Bound +
Footnote †: journal: Sensors
## 4 The DCPVP Protocol
### Phase B: Cluster Head Election Phase  
In this phase, each node calculates its own weight:  
\[W\ =\ \alpha\ F(\text{MD})\ +\ \beta\ F(\Delta C)\ +\ \gamma\ F(Y)\ -\ \phi\ F(\frac{E_{res}}{E_{0}}),0\leq\alpha,\beta,\gamma,\phi\ \leq 1 \tag{8}\]  
where \(\alpha\), \(\beta\), \(\gamma\) and \(\phi\) are adjust coefficients, Y is the times that the node has been a CH so far, \(E_{res}\) is the residual energy, \(E_{0}\) is initial energy, MD is mean distance to neighbors and \(\Delta C\) is the optimum deviation. MD can be calculated as follows:  
\[MD\ =\frac{\sum\limits_{i=1}^{N_{S}}\text{DV(i,x)}}{N_{s}} \tag{9}\]  
where X is the nodes' ID, \(N_{s}\) is the number of neighbors, DV is the distance vector, DV (i, j) is the distance between the nodes i and j and \(\Delta C\) can be calculated as follows:  
\[\Delta C\ =\ |\ N_{s}\ -\ N_{0}\ | \tag{10}\]  
where \(N_{s}\) is the number of neighbors and \(N_{O}\) is optimum number of the neighbors (See Figure 2):  
\[N_{O}\ =\ floor\ (N_{r}) \tag{11}\]  
\[N_{r}\ =\ \left\{\begin{array}{ll}1\ \times\ Nm&D>Dt1\\ C1\ \times Nm&Dt1>D>Dt2\\ C2\ \times Nm&Dt2>D>Dt3\\ C3\ \times Nm&Dt3>D\end{array}\right.,\quad 0<C3<C2<C1<1 \tag{12}\]  
where \(D\) is the distance from BS, \(Nm\) is maximum cluster size, \(Dt1\), \(Dt2\) and \(Dt3\) are threshold values of distance and \(C1\), \(C2\) and \(C3\) are the coefficients and are less than 1 (\(C1\), \(C2\),\(C3<1\)). Thus we allow the nodes which are farther from BS to form bigger clusters, and allow the closer nodes to only form smaller clusters. This allows load balancing. In this step, each node calculates its own weight and broadcasts it sequentially. Now each node adds a column to its table and writes the weight of each neighbor in it as Table 2.  
Then each node creates a priority list from its neighbors based on their weight. After that they broadcast a packet containing the voting list (Vote-Packet) and vote the best nodes. In other words, the Vote-Packet is a list of nodes' IDs which are sorted based on their weight. After finishing the election, the node which has the most number of votes, is chosen as the CH and introduces itself by sending a packet (CH-Packet). Figures 5 and 6 show the flowchart of this phase.  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Neighbors** & **ID** & **Distance (m)** & **Weight** \\ \hline
1 & 32 & 0.56 & W\({}_{1}\) \\
2 & 57 & 4.41 & W\({}_{2}\) \\... &... &... &... \\ M & 4 & 3.1 & W\({}_{m}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Nodes information from network topology in every round.  
All ID-based sequential broadcasts have the same following algorithm. When node i wants to send its own packet, it waits for i-1 time slots from the end of last step until its turn comes and then sends its own packet. During all other time slots the node monitors the channel for receiving probable packets [33].
# A Dual Stage Low Power Converter Driving for Piezoelectric Actuator Applied in Micro Mobile Robot
## 3 Design and Control of the Drive Circuit  
A dual-stage circuit is designed to meet the following requirements: The first stage should step up the low battery voltage to a high voltage, whereas the second stage should transform the DC voltage into a time-varying drive signal.  
The drive circuit topology of the piezoelectric bimorph proposed in this study is depicted in Figure 4. It includes a flyback converter as the DC/DC stage and a bidirectional active half-bridge converter as the high-voltage DC/AC drive stage [28]. The flyback converter boosts the low voltage of the lithium battery to the DC high voltage required for the driving actuator. The power metal-oxide-semiconductor field-effect transistor (MOSFET) in both stages is controlled, via pulse width modulation (PWM), in which a novel control strategy is utilized in a half-bridge drive stage to output unipolar signals of any shape to satisfy the actuator drive requirements. The selected inductor with a high quality factor in the half-bridge can effectively achieve energy transfer. The half-bridge structure has less footprint area than other DC/AC converters [29], thereby increasing the energy density of the total system.  
Figure 3: Driving methods for piezoelectric actuators: (**a**) Bimorph with alternating driving; (**b**) bimorph with simultaneous driving.
# Solar Radiation Pressure Models for BeiDou-3 I2-S Satellite: Comparison and Augmentation
###### Abstract  
As one of the most essential modeling aspects for precise orbit determination, solar radiation pressure (SRP) is the largest non-gravitational force acting on a navigation satellite. This study focuses on SRP modeling of the BeiDou-3 experimental satellite I2-S (PRN C32), for which an obvious modeling deficiency that is related to SRP was formerly identified. The satellite laser ranging (SLR) validation demonstrated that the orbit of BeiDou-3 I2-S determined with empirical 5-parameter Extended CODE (Center for Orbit Determination in Europe) Orbit Model (ECOM1) has the sun elongation angle (\(\varepsilon\) angle) dependent systematic error, as well as a bias of approximately \(-16.9\) cm. Similar performance has been identified for European Galileo and Japanese QZSS Michibiki satellite as well, and can be reduced with the extended ECOM model (ECOM2), or by using the a priori SRP model to augment ECOM1. In this study, the performances of the widely used SRP models for GNSS (Global Navigation Satellite System) satellites, i.e., ECOM1, ECOM2, and adjustable box-wing model have been compared and analyzed for BeiDou-3 I2-S satellite. In addition, the a priori SRP models are derived based on analytical cuboid box model and empirically spectra analysis, respectively. Use of the a priori model combined with ECOM1 was finally demonstrated to reduce the \(\varepsilon\)-angle-dependent systematic error, and thus improved the radial orbit accuracy by nearly 35 per cent when compared to the solution with standalone ECOM1, as revealed by the one way SLR residuals.  
BeiDou-3 I2-S; precise orbit determination; solar radiation pressure; box-wing Article
# Modeling GPR Wave Propagation in Complex Underground Structures Using Conformal ADI-FDTD Algorithm
## 2 Methodology
### UPML Absorbing Boundary Condition
 \tag{44}\]  
In sub-step \(n{\rightarrow}n+1/2\), discretize equation (38) as  
\[D_{z}^{n+\frac{1}{2}}(i,j)=\frac{bx}{ax}.D_{z}^{n}(i,j)+\frac{1}{ax}.\bigg{\{} \frac{1}{ax}.\bigg{[}H_{y}^{n+\frac{1}{2}}\bigg{(}i+\frac{1}{2},j\bigg{)}-H_{ y}^{n+\frac{1}{2}}\bigg{(}i-\frac{1}{2},j\bigg{)}\bigg{]}-\frac{1}{ \Delta y}.\bigg{[}H_{x}^{n}\bigg{(}i,j+\frac{1}{2}\bigg{)}-H_{x}^{n}\bigg{(}i,j-\frac{1}{2}\bigg{)}\bigg{]}\bigg{\}}. \tag{45}\]  
The discretization scheme of Equation (39) is  
\[E_{z}^{n+\frac{1}{2}}(i,j)=\frac{by}{ay}.E_{z}^{n}(i,j)+\frac{er}{ay}.\bigg{[} D_{z}^{n+\frac{1}{2}}(i,j)-D_{z}^{n}(i,j)\bigg{]}, \tag{46}\]  
where  
\[D_{z}^{n+\frac{1}{2}}(i,j)=D_{z}^{n}(i,j)+\frac{ay}{er}.E_{z}^{n+\frac{1}{2}}( i,j)-\frac{by}{er}.E_{z}^{n}(i,j). \tag{47}\]  
The discretization schemes of Equations (40)-(43), respectively, are  
\[B_{x}^{n+\frac{1}{2}}\bigg{(}i,j+\frac{1}{2}\bigg{)}=\frac{by}{ay}.B_{x}^{n} \bigg{(}i,j+\frac{1}{2}\bigg{)}-\frac{1}{ay}.\frac{1}{\Delta y}.[E_{z}^{n}(i,j +1)-E_{z}^{n}(i,j)], \tag{48}\]  
\[B_{y}^{n+\frac{1}{2}}\bigg{(}i+\frac{1}{2},j\bigg{)}=\frac{bx}{ax}.B_{y}^{n} \bigg{(}i+\frac{1}{2},j\bigg{)}+\frac{1}{ax}.\frac{1}{\Delta x}.\bigg{[}E_{z}^ {n+\frac{1}{2}}(i+1,j)-E_{z}^{n+\frac{1}{2}}(i,j)\bigg{]}, \tag{49}\]  
\[H_{x}^{n+\frac{1}{2}}\bigg{(}i,j+\frac{1}{2}\bigg{)}=H_{x}^{n}\bigg{(}i,j+ \frac{1}{2}\bigg{)}+hr.\bigg{[}ax.B_{x}^{n+\frac{1}{2}}\bigg{(}i,j+\frac{1}{2} \bigg{)}-bx.B_{x}^{n}\bigg{(}i,j+\frac{1}{2}\bigg{)}\bigg{]}, \tag{50}\]  
\[H_{y}^{n+\frac{1}{2}}\bigg{(}i+\frac{1}{2},j\bigg{)}=H_{y}^{n}\bigg{(}i+\frac{1 }{2},j\bigg{)}+hr.\bigg{[}ay.B_{y}^{n+\frac{1}{2}}\bigg{(}i+\frac{1}{2},j\bigg{)} \bigg{)}-by.B_{y}^{n}\bigg{(}i+\frac{1}{2},j\bigg{)}\bigg{]}, \tag{51}\]  
In Equations (45)-(51), except for Equation (48), both sides contain the field component at time (\(n+1/2\)). By eliminating the simultaneous components from Equations (45)-(51), except for \(E_{z}{}^{n+1/2}\), we obtain  
\[a_{i}E_{z}^{n+\frac{1}{2}}(i-1,j)+b_{i}E_{z}^{n+\frac{1}{2}}(i,j)+c_{i}E_{z}^{ n+\frac{1}{2}}(i+1,j)=d_{i}, \tag{52}\]  
where\[ \begin{array}{l}a_{i}=-\frac{hr}{(\Delta x)^{2}}\cdot\frac{a\chi}{ax}\quad b_{i}= \frac{ax\cdot xy}{er}+\frac{hr}{(\Delta x)^{2}}\cdot\frac{2a\chi}{ax}\quad c_{i}= -\frac{hr}{(\Delta x)^{2}}\cdot\left[ay\Big{(}i+\frac{1}{2},i\Big{)}\right/ax \Big{(}i+\frac{1}{2},j\Big{)}\right]\\ d_{i}=\frac{ax\cdot yy}{er}E_{z}^{n}(i,j)+\Big{\{}\frac{1}{\Delta x}\cdot\left[H_ {y}^{n}\Big{(}i+\frac{1}{2},j\Big{)}-H_{y}^{n}\Big{(}i-\frac{1}{2},j\Big{)} \right]-\frac{1}{\Delta y}\cdot\left[H_{x}^{n}\Big{(}i,j+\frac{1}{2}\Big{)}-H _{x}^{n}\Big{(}i,j-\frac{1}{2}\Big{)}\right]\Big{\}}\\ +(bx-ax)\cdot D_{z}^{n}(i,j)+\frac{hr}{\Delta x}\cdot\left[\frac{ay\cdot bx}{ ax}-by\right]\cdot B_{y}^{n}\Big{(}i+\frac{1}{2},j\Big{)}-\left[\frac{ay\cdot bx}{ ax}-by\right]\cdot B_{y}^{n}\Big{(}i-\frac{1}{2},j\Big{)}.\end{array}