# File Overview and Run Order

This repository contains generated data and plotting scripts for the paper
figures. Most plotting scripts can be run directly if the intermediate data
already exist. If data are missing, run the corresponding data-generation
scripts first.

Recommended working rule:

```bash
cd <project_root>
conda activate twistpy_2024_tf
python <script.py>
```

Use the `seismology` environment for scripts that generate synthetic waveforms
with `pyGRT`. Use `twistpy_2024_tf` for polarization analysis, data processing,
and plotting unless a script note says otherwise.


# Required Environments

Two Anaconda environments are used:

- `seismology`: includes `pyGRT`; used for synthetic waveform generation. The requirements shown in `seismology_requirements.yml`.
- `twistpy_2024_tf`: used for polarization analysis, waveform processing, FFT,
  STFT, pressure/wind analysis, and plotting.  The requirements shown in `twistpy_2024_tf_requirements.yml`


Important shared files:

- `twistpy_functions.py` and `foundations.py`: utility functions used by the
  plotting and analysis scripts.
- `MQget.py`: script for downloading InSight marsquake data from IRIS.
    If `./InsightMarsquakeData/MarsData_Sol/` not existed, use Part A in script to download.
    If `./InsightMarsquakeData/Marsdata_Deck/` not existed, use Part C.
    If `./InsightMarsquakeData/MarsquakeLonger/` not existed, use Part B.
- `XB.ELYSE.BH.response.xml`, `XB.ELYSE.BH.response.all_20250701.xml`, and
  `SP_XB.ELYSE.HH.response_all.xml`: instrumental response files.
- `events_mars_extended_multiorigin_v14_2023-01-01.xml`: Marsquake catalog.


!!!!!!!
Take Care:
The data have been generated. Don't need to run scripts with '_datagenerate' to save time.

If you want to generate by yourself, run scripts in `Global Pre-Run Checklist` first.

# Global Pre-Run Checklist
Before generating the final paper figures, make sure the following intermediate
products exist or are generated in this order:

1. All-Sol polarization tables:
Take care: The generation will cost a lot of time. 

```bash
cd Figure_with_BAZ_INC_EFFR_statics
python ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate.py
```
or use `ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate_multi.py` for
multi processing at the same time. Default use 20 kernel cpu.


Outputs:

- `Figure_with_BAZ_INC_EFFR_statics/Figures/BAZ_table.csv`
- `Figure_with_BAZ_INC_EFFR_statics/Figures/INC_table.csv`
- `Figure_with_BAZ_INC_EFFR_statics/Figures/EFFR_table.csv`

Table S2 and Table S3 was from these files.

These tables are used by `Main_Figure_plot/Fig2_2_4_feature.py`,
`Main_Figure_plot/Fig3_forward_for_BAZ_INC_20260904.py`,
`Figure_with_Synthetics/FigS22_force_changed_BAZ_S2_slice.py`.

2. Event FFT peak table:

```bash
cd Figure_with_event_FFT_STFT
python ExtendedDataFig6_and_FigS23_event_FFT_powerRatio.py
```

Output used by Fig. 1:

- `Figure_with_event_FFT_STFT/2_4_PR_envelop_withlandermodes/P120s_fft_envelope_peak_values.csv`

Note: the script currently has `plot_scatter_from_csv_only = True`. If the CSV
does not exist, set it to `False` or use the existing CSV.

3. Wind-direction KDE tables:

```bash
cd Pressure
python Fig3_2_4_BAZ_wind_direction_kde_BAZ_datagenerate.py
```

Outputs:

- `Pressure/wind_direction_kde_tables/BMY_WIND_DIRECTION_KDE_1730_2230_table.csv`
- `Pressure/wind_direction_kde_tables/BPY_WIND_DIRECTION_KDE_1730_2230_table.csv`

These tables are required before running
`Main_Figure_plot/Fig3_forward_for_BAZ_INC_20260904.py`.

4. Force-changed synthetic model products for Fig. 3:

```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig3_1_Synthetic_Models_statistics_force_changed_datagenerate.py

conda activate twistpy_2024_tf
python Fig3_2_Synthetic_Models_statistics_force_changed_polaranalysis_generate.py
python Fig3_3_force_changed_peak_4view_save_csv.py
```

Outputs:

- `Figure_with_Synthetics/Lander_Force_Test/Synthetic_2_4Hz_Model_Hobiger2/ratio_*/2htest_clean.mseed`
- `Figure_with_Synthetics/Lander_Force_Test/Synthetic_2_4Hz_Model_Hobiger2/ratio_*/2htest_noisy.mseed`
- `Figure_with_Synthetics/Lander_Force_Test/Figures/ExtendDataFigX_force_changed_polarization_peak_sigma_window2.5s_101019.csv`

5. Other comparison-model synthetic products for Fig. 3:

```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig3_4_Synthetic_Models_statistics_for_Xiao_Hobiger_BMCBMN_datagenerate.py
```

This script should be run for the required comparison batches, for example the
default model folders (actually`_101010` but have no suffix) 
and the `_101019` model folders used by
`Main_Figure_plot/Fig3_forward_for_BAZ_INC_20260904.py`.

6. Vp/Vs grid products for Fig. 4:

```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig4_de_VelocityModels_statistics_datagenerate1_modelsynthetic.py

conda activate twistpy_2024_tf
python Fig4_de_VelocityModels_statistics_datagenerate2.py
```

Important check: `Fig4_de_VelocityModels_statistics_datagenerate1_modelsynthetic.py`
For the main Fig. 4 workflow, make sure the
generated folders and files are in `Figure_with_Synthetics/Layer_Model_VpVs_forward_20260903/`
or update the reader paths consistently.


# Main Figure Scripts

Run the final main-figure scripts from `Main_Figure_plot` after the global
pre-run checklist is complete.

## Fig. 1

```bash
conda activate twistpy_2024_tf
cd Main_Figure_plot
python Fig1_why_24_important.py
```

Required before running:
- `Main_Figure_plot/Figure/Figure1_hor_20260909.png`
- `InsightMarsquakeData/MarsData_Sol/`
- `Figure_with_event_FFT_STFT/2_4_PR_envelop_withlandermodes/P120s_fft_envelope_peak_values.csv`

Output:
- `Main_Figure_plot/Figure/Fig1_why_24_import.png`

## Fig. 2

```bash
conda activate twistpy_2024_tf
cd Main_Figure_plot
python Fig2_2_4_feature.py
```

Required before running:
- `Figure_with_BAZ_INC_EFFR_statics/Figures/BAZ_table.csv`
- `Figure_with_BAZ_INC_EFFR_statics/Figures/INC_table.csv`
- required Sol 793 polarization `.npy` and waveform files in
  `Figure_with_BAZ_INC_EFFR_statics/BAZ_data/`

If these are missing, first run:
```bash
cd Figure_with_BAZ_INC_EFFR_statics
python ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate.py
```
or `ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate_multi.py`

Output:

- `Main_Figure_plot/Figure/Fig2_24_feature.png`

## Fig. 3 and Table 1 & 2 (Xiao2022, Hobiger1, Hobiger2 results)

Run:
```bash
conda activate twistpy_2024_tf
cd ../Main_Figure_plot
python Fig3_forward_for_BAZ_INC_20260904.py
```

If something is empty, refer Global Pre-Run Checklist 4 and 5, and run 
```bash
cd ../Pressure
conda activate twistpy_2024_tf
python Fig3_2_4_BAZ_wind_direction_kde_BAZ_datagenerate.py
```
To generate Fig3 Sol 301 wind direction KDE result.

Required before final plotting:
- `Figure_with_Synthetics/Fig4_Models_compare/Synthetic_2_4Hz_Model_*/`
- `Figure_with_Synthetics/Fig4_Models_compare/Synthetic_2_4Hz_Model_*_101019/`
- `Figure_with_Synthetics/Lander_Force_Test/Figures/ExtendDataFigX_force_changed_polarization_peak_sigma_window2.5s.csv`
- `Figure_with_BAZ_INC_EFFR_statics/Figures/BAZ_table.csv`
- `Pressure/wind_direction_kde_tables/BMY_WIND_DIRECTION_KDE_1730_2230_table.csv`
- `Pressure/wind_direction_kde_tables/BPY_WIND_DIRECTION_KDE_1730_2230_table.csv`

Output:

- `Main_Figure_plot/Figure/Fig3_forward_BAZ_INC_Hobiger2_noisy_BAZ_slice.png`

## Fig. 4 and Table 2 (BMC, BMN results)

Run:
```bash
conda activate twistpy_2024_tf
cd ../Main_Figure_plot
python Fig4_inversion_for_subsurface_vpvs.py
```

If something empty, run:
```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig4_de_VelocityModels_statistics_datagenerate1_modelsynthetic.py

conda activate twistpy_2024_tf
python Fig4_de_VelocityModels_statistics_datagenerate2.py
```

Required before final plotting:

- `Figure_with_Synthetics/Layer_Model_VpVs_forward_20260903/Synthetic_2_4Hz_Model_*/2htest_clean.mseed`
- `Figure_with_Synthetics/Layer_Model_VpVs_forward_20260903/Synthetic_2_4Hz_Model_*/2htest_noisy.mseed`
- corresponding `cleandata/` and `noisedata/` polarization `.npy` files

Output:
- `Main_Figure_plot/Figure/Fig4_inversion_for_subsurface_vpvs_2.5s_ratio_1p0_1p0_1p9.png`


# Figure_with_BAZ_INC_EFFR_statics

Directory contents:

- `BAZ_data/`: polarization-analysis `.npy` and waveform products.
- `Figures/`: Extended Data Fig. 3, Extended Data Fig. 7, and CSV tables.
- `Figures_S10-15/`: Supplementary Fig. 10-15 products if generated.

Run order:

1. Generate all-Sol BAZ/INC/EFFR tables and required `.npy` files:
Same as `Global Pre-Run Checklist 1`.
```bash
# no need to run if data existed. 
cd Figure_with_BAZ_INC_EFFR_statics
python ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate.py 
```

2. Plot Extended Data Fig. 3:

```bash
python ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR.py
```

3. Generate window-length test products for Extended Data Fig. 7:

```bash
# no need to run if data existed.
python ExtendDataFig7_window_test_datagenerate.py  

# Plot figure.
python ExtendDataFig7_window_test.py
```

4. Plot Supplementary Fig. 10-15:

```bash
python FigS10_15_Landermodes_BAZ_INC_EFFR_statics.py
```

Note: `FigS10_15_Landermodes_BAZ_INC_EFFR_statics.py` reads files from
`BAZ_data/`. If the needed Sol/frequency/window `.npy` files are missing, run
the all-Sol data-generation script first or set `BAZcompute=True` in the script.


# Figure_with_event_FFT_STFT

Directory contents:

- `FFT/`: FFT results for impact events.
- `Marsquake_before_after_data/`: waveform windows before and after P arrival.
- `STFT/`: STFT figures.
- `2_4_PR_envelop_withlandermodes/`: FFT-envelope peak figures and CSV table with lander modes.
- `2_4_PR_envelop/`: FFT-envelope peak figures and its CSV table without lander modes, refer to Text S2.

Run order:

1. Generate or plot the FFT-envelope peak products:
Before run, select corresponding parameters during line number 384 - 416 in script. 
```bash
cd Figure_with_event_FFT_STFT
python ExtendedDataFig6_and_FigS23_event_FFT_powerRatio.py
```
Output:

- `2_4_PR_envelop_withlandermodes/P120s_fft_envelope_peak_values.csv`
- Extended Data Fig. 6 / Fig. S23 FFT-envelope scatter plots.

2. Plot FFT before-after comparisons for nine impact events:

```bash
python FigS21_fft_after_minus_before_9events_multiband.py
```

3. Plot event/Sol STFT figures:

```bash
python FigS_STFT_of_Event_and_Sol_20260507.py
python FigS1_2_STFT_of_S1000a_S1094b.py
python FigS1_2_STFT_of_other_events.py
```

# Figure_with_polar_motion

Directory contents:

- `Figures/`: figures created by the scripts below.
- `mars_relief_cache/`: GMT Mars map cache.
- `Polar_motion_data/`: filtered and sliced waveform data for impact events.
- `Figures_Polarcheck_longer/`: filtered and sliced waveform and polarization 
    motion figures for VF impact events -10s to 10s from P arrival.
- `Figures_Polarcheck/`: filtered and sliced waveform and polarization 
    motion figures for VF impact events for selected time range.
- `mars_colorbar.cpt`: color map used by the Mars map script.

Run order:

1. Generate Mars map data:

```bash
cd Figure_with_polar_motion
conda activate seismology
python ExtendDataFigure1_MarsMap_datagenerate.py
```

2. Plot S1000a/S1094b polar motion:

```bash
conda activate twistpy_2024_tf
python ExtendDataFig2_S1000aS1094b_polarmotion.py
```

3. Plot impact-event filter-band polar motion figures:

```bash
python FigS3_9_ImpactAllFilterBand.py
```

Optional/deprecated:
This script plot -10 to 10s from P arrival for every impact event. 
- `FigS3_9_ImpactAllFilterBand_longer.py`


# Figure_with_Synthetics

Directory contents:

- `BMN_Simulated_waveform_landermodes/`: lander modes synthetic products.
- `Fig4_Models_compare/`: synthetic waveforms and polarization products for
  Hobiger1, Hobiger2, and Xiao2022 models. Actually for Fig. 3
- `Fig4_Models_compare_testrandomforce/` make sure the random seed not influence the  
    results. Actually for Fig. 3
- `Lander_Force_Test/`: force-ratio synthetic tests for Fig. 3.
- `Layer_Model_VpVs_forward_20260903/`: Vp/Vs perturbation synthetic tests for Fig. 4.
- `Figures/`: synthetic-model figures and summary CSV tables.

## Fig. 3 Synthetic Workflow
Refer to Global Pre-Run Checklist 4 and 5.

If data not exist, run in this order:
```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig3_1_Synthetic_Models_statistics_force_changed_datagenerate.py
conda activate twistpy_2024_tf
python Fig3_2_Synthetic_Models_statistics_force_changed_polaranalysis_generate.py
python Fig3_3_force_changed_peak_4view_save_csv.py
conda activate seismology
python Fig3_4_Synthetic_Models_statistics_for_Xiao_Hobiger_BMCBMN_datagenerate.py
```

If data existed,
Then run the main Fig. 3 script:
```bash
cd ../Main_Figure_plot
conda activate twistpy_2024_tf
python Fig3_forward_for_BAZ_INC_20260904.py
```

Additional Fig. 3 related supplementary plot:

```bash
cd ../Figure_with_Synthetics
python FigS22_force_changed_BAZ_S2_slice.py
```

`FigS22_force_changed_BAZ_S2_slice.py` requires:

- `Lander_Force_Test/Figures/ExtendDataFigX_force_changed_polarization_peak_sigma_window2.5s_101019.csv`
- `Figure_with_BAZ_INC_EFFR_statics/Figures/BAZ_table.csv`

## Fig. 4 Synthetic Workflow
Refer to Global Pre-Run Checklist 6.

If data did not existed.
Run in this order:
```bash
cd Figure_with_Synthetics
conda activate seismology
python Fig4_de_VelocityModels_statistics_datagenerate1_modelsynthetic.py

conda activate twistpy_2024_tf
python Fig4_de_VelocityModels_statistics_datagenerate2.py
```

Then run:

```bash
cd ../Main_Figure_plot
python Fig4_inversion_for_subsurface_vpvs.py
```

- `Figure_with_Synthetics/Figures/Fig4_Model_compare_polarization_peak_sigma_window2.5s_VpVs_grid_ratio_1p0_1p0_1p9.csv`

## Lander-Mode Synthetic Workflow

Run in this order:

```bash
cd Figure_with_Synthetics
conda activate seismology
# If data existed, skip this.
python ExtendDataFig5_FigS17_20_Synthetic_Models_statistics_for_LanderModes_datagenerate.py


conda activate twistpy_2024_tf
python ExtendDataFig5_FigS17_20_Synthetic_Models_statistics_for_LanderModes.py
```

Important check: the data-generation script currently assigns `savefigdir` to
`./testdir/` after assigning `./BMN_Simulated_waveform_landermodes/`. If you want
the plotting script to read `BMN_Simulated_waveform_landermodes/`, make sure the
generated `results/` files are there or update both scripts consistently.

## Random-Force Test Workflow
Used to comfirm the numpy random force have no influence to the results.
Run in this order:
```bash
cd Figure_with_Synthetics
conda activate seismology
python EDFigSX_Synthetic_Models_statistics_test_randomforce_datagenerate.py
conda activate twistpy_2024_tf
python EDFigSX_Synthetic_Models_statistics_test_randomforce_plot.py
```

Outputs are written under:
- `Figure_with_Synthetics/Fig4_Models_compare_testrandomforce/`


# Pressure

Directory contents:

- `InSight_PS_TWINS_csv/`: pressure and TWINS wind CSV files.
- `Figures/`: pressure, wind, H/V, and BAZ figures.
- `Figures/HVwind/`: Extended Data Fig. 4 H/V-wind-BAZ products.
- `wind_direction_kde_tables/`: wind-direction KDE peak/sigma tables.
- `ps_twins_data_download_select_complete_data_file.py ` and `twins_data_download_all.py/`: 
    download all TWINS and Pressure data.
    `InSight_PS_TWINS_csv` did not contain all data because the storage. If nee to use other data, use these two scripts to download or use:
    https://atmos.nmsu.edu/PDS/data/PDS4/InSight/twins_bundle/data_calibrated/  for TWINS
    and 
    https://atmos.nmsu.edu/PDS/data/PDS4/InSight/ps_bundle/data_calibrated/ for Pressure.


Run order:

1. Generate wind-direction KDE tables needed by Fig. 3:

```bash
cd Pressure
python Fig3_2_4_BAZ_wind_direction_kde_BAZ_datagenerate.py
```

Just use functions from these two scripts.
  Plot Sol 298 pressure/wind/HV/STFT ratio:

  ```bash
  python plot_sol298_pressure_wind_hv_stft_ratio.py
  ```
  Plot Sol 298 wind and BAZ LMST panels:
  ```bash
  python plot_sol298_wind_baz_lmst.py
  ```


2. Plot Extended Data Fig. 4 H/V-wind-BAZ statistics and Table S4:
For Sol 298 and 301(change Target_Sol parameter in script).
```bash
python ExtendDataFig4_LMST_Sol298_wind_baz_statistic.py
```

`plot_sol298_wind_baz_lmst.py` and `ExtendDataFig4_LMST_Sol298_wind_baz_statistic.py`
read BAZ/INC/EFFR `.npy` products from

`Figure_with_BAZ_INC_EFFR_statics/BAZ_data/`. Generate them first with:
```bash
cd ../Figure_with_BAZ_INC_EFFR_statics
python ExtendDataFig3_mainpaper_allSolsPlotBAZ_INC_EFFR_datagenerate.py
```


# InsightMarsquakeData

Directory contents:

- `deglitch_bash_sol/`: Seisglitch scripts to deglitch data, use SeisGlitch python package.
- `ELYS0_68_SH/`: EH and SH channel files.
- `MarsData_Sol/`: InSight Mars data organized by Sol.
- `MarsquakeLonger/`: marsquake data for the nine InSight impact marsquakes.

No figure should start from this directory directly. It is used as the waveform
input directory by the plotting and data-generation scripts above.

- `impact_event_time.py` output the impact events P and S arrival time.

# UsedPackages

Directory contents:

- `insight-seis-catalog/`
- `marstimeconverter/`

These local packages are added to `sys.path` by most scripts. If imports fail,
install them from source or keep the directory layout unchanged.


# Some of Reference Python Packages
These packages in `./UsedPackages/`

## Twistpy
We use this package for polarization filter.
https://twistpy.org/installation.html

## marstimeconverter
https://gitlab.obspm.fr/gsainton/marstimeconverter/-/wikis/Installation-and-Setup

Better to Installation from source (for development).

We need to use the development folder to get config.
When we use python script, we need to use these codes:

```python
from marstimeconverter.converter import MarsTimeConverter
from marstimeconverter.plotting import plot_lmst_vs_utc
config = "./marstimeconverter/CONFIG/mission_config_file.xml"
mtc = MarsTimeConverter(mission_config_file=config, mission="InSight")
```

## insight-seis-catalog
We use this package to read the insight catalog.
https://github.com/sceylan/insight-seis-catalog/tree/master

## SeisGlitch
We use this package to deglitch insight marsquake data.
https://seisglitch.readthedocs.io/en/latest/