# PySME how to 

This page describes some new (and in testing) function of PySME.

## How to get the atmosphere grid

```py
from pysme.sme import SME_Structure
from pysme.synthesize import Synthesizer

# Set up sme structure
sme = SME_Structure()
# Define stellar parameters
sme.teff, sme.logg, sme.monh, sme.vmic = 5777, 4.485, 0, 1.46
# Select atmosphere source; optional, default is MARCS models ('marcs2014.sav')
# Available options are listed in https://pysme-astro.readthedocs.io/en/latest/concepts/lfs.html
sme.atmo.source = 'marcs2014.sav'
synthesizer = Synthesizer()
# Interplate the atmopshere
atmo = synthesizer.get_atmosphere(sme).atmo
```

The atmosphere grid is stored in `atmo`:

```py
Atmosphere(teff=5777.0,
logg=4.485,
abund= [M/H]=0.000 applied to abundance pattern. Values below are abundances.
  H      He     Li     Be     B      C      N      O      F      Ne     Na   
 12.000 12.114  1.046  1.376  2.696  8.386  7.776  8.656  4.556  7.836  6.166
  Mg     Al     Si     P      S      Cl     Ar     K      Ca     Sc     Ti   
  7.526  6.366  7.506  5.356  7.136  5.496  6.176  5.076  6.306  3.166  4.896
  V      Cr     Mn     Fe     Co     Ni     Cu     Zn     Ga     Ge     As   
  3.996  5.636  5.386  7.446  4.916  6.226  4.206  4.596  2.876  3.576  2.286
  Se     Br     Kr     Rb     Sr     Y      Zr     Nb     Mo     Tc     Ru   
  3.326  2.556  3.246  2.596  2.916  2.206  2.576  1.416  1.916 -7.964  1.836
  Rh     Pd     Ag     Cd     In     Sn     Sb     Te     I      Xe     Cs   
  1.116  1.656  0.936  1.766  1.596  1.996  0.996  2.186  1.506  2.236  1.066
  Ba     La     Ce     Pr     Nd     Pm     Sm     Eu     Gd     Tb     Dy   
  2.166  1.126  1.696  0.576  1.446 -7.964  0.996  0.516  1.106  0.276  1.136
  Ho     Er     Tm     Yb     Lu     Hf     Ta     W      Re     Os     Ir   
  0.506  0.926 -0.004  1.076  0.056  0.876 -0.174  1.106  0.226  1.246  1.376
  Pt     Au     Hg     Tl     Pb     Bi     Po     At     Rn     Fr     Ra   
  1.636  1.006  1.126  0.896  1.996  0.646 -7.964 -7.964 -7.964 -7.964 -7.964
  Ac     Th     Pa     U      Np     Pu     Am     Cm     Bk     Cf     Es   
 -7.964  0.056 -7.964 -0.524 -7.964 -7.964 -7.964 -7.964 -7.964 -7.964 -7.964,
vturb=1.0,
lonh=1.5,
source='marcs2012.sav',
method='grid',
geom='PP',
radius=1.0,
height=array([-3.28991027e+10,  1.08125362e+01,  7.27886834e+16, ...,
       -3.07627183e+21,  1.08125000e+01,  6.88842390e+27]),
opflag=array([1, 1, 1, ..., 0, 0, 0]),
wlstd=5000.0,
depth='RHOX',
interp='TAU',
rhox=array([0.01206928, 0.0155857 , 0.02009963, ..., 6.73889351, 7.1708134 ,
       7.68721481]),
tau=array([1.89675940e-05, 2.95631869e-05, 4.58225743e-05, ...,
       2.44857945e+01, 3.79713691e+01, 5.92142993e+01]),
temp=array([4101.87925384, 4145.25143759, 4189.80397932, ..., 9422.35229943,
       9677.82956675, 9935.68732234]),
rho=array([1.36430747e-09, 1.74484534e-09, 2.22803405e-09, ...,
       3.24239421e-07, 3.33732002e-07, 3.45981912e-07]),
xna=array([6.48271628e+14, 8.29097296e+14, 1.05869036e+15, ...,
       1.54084995e+17, 1.58596875e+17, 1.64417830e+17]),
xne=array([5.05875737e+10, 6.47950015e+10, 8.29144374e+10, ...,
       4.09172378e+15, 5.26581938e+15, 6.73330283e+15]),
citation_info='\n                @ARTICLE{2008A&A...486..951G,\n                    author = {{Gustafsson}, B. and {Edvardsson}, B. and {Eriksson}, K. and\n                    {J{\\o}rgensen}, U.~G. and {Nordlund}, {\\r{A}}. and {Plez}, B.},\n                    title = "{A grid of MARCS model atmospheres for late-type stars. I. Methods and general properties}",\n                    journal = {Astronomy and Astrophysics},\n                    keywords = {stars: atmospheres, Sun: abundances, stars: fundamental parameters, stars: general, stars: late-type, stars: supergiants, Astrophysics},\n                    year = "2008",\n                    month = "Aug",\n                    volume = {486},\n                    number = {3},\n                    pages = {951-970},\n                    doi = {10.1051/0004-6361:200809724},\n                    archivePrefix = {arXiv},\n                    eprint = {0805.0554},\n                    primaryClass = {astro-ph},\n                    adsurl = {https://ui.adsabs.harvard.edu/abs/2008A&A...486..951G},\n                    adsnote = {Provided by the SAO/NASA Astrophysics Data System}\n                }\n            ')
```

## How to update central depth and line range

Two new parameters for each line, central depth and line range, are added to PySME since 0.4.189. 

Central depth is the line depth **before** any kind of boradening.
It should **not** be used as an indicator of how deep a line is in the final synthetic spectra, since actual line depth also depands on the broadning mechanisms.
However, it can be used for screnning out the weak lines in the line list, thus acclerating the synthesis speed.

Line range indicates the wavelength range where a line has contribution, and it can be used to select the revelant lines for the synthetic wavelength, also for acclerating the synthesis speed.

These parameters can be calculated as follow:

```py
from pysme.sme import SME_Structure
from pysme.linelist.vald import ValdFile
from pysme.synthesize import synthesize_spectrum, Synthesizer
from pysme.solve import solve

linelist = ValdFile('SME/examples/sun.lin')
```

After we read in the line list, it looks like:

```py
   species     wlcent  gflog     excit  j_lo     e_upp  j_up  lande_lower  \
0     Fe 1  6436.4055 -2.460  4.186364   2.0  6.112128   1.0         0.68   
1     Eu 2  6437.6106 -1.242  1.319612   5.0  3.245008   5.0         1.73   
2     Eu 2  6437.6121 -1.280  1.319636   5.0  3.245027   5.0         1.73   
3     Eu 2  6437.6135 -2.473  1.319612   5.0  3.245008   5.0         1.73   
4     Eu 2  6437.6194 -2.511  1.319636   5.0  3.245027   5.0         1.73   
..     ...        ...    ...       ...   ...       ...   ...          ...   
62    Mn 1  6443.4689 -2.818  3.772307   1.5  5.695960   2.5         1.21   
63    Mn 1  6443.4726 -3.538  3.772307   1.5  5.695960   2.5         1.21   
64    Mn 1  6443.4821 -2.275  3.772307   1.5  5.695960   2.5         1.21   
65    Mn 1  6443.4887 -2.964  3.772307   1.5  5.695960   2.5         1.21   
66    Mn 1  6443.4938 -3.918  3.772307   1.5  5.695960   2.5         1.21   

    lande_upper  lande  ...  gamvw  depth  \
0          0.56   0.73  ... -7.810  0.112   
1          1.79   1.76  ...  0.000  0.011   
2          1.79   1.76  ...  0.000  0.010   
3          1.79   1.76  ...  0.000  0.011   
4          1.79   1.76  ...  0.000  0.010   
..          ...    ...  ...    ...    ...   
62         1.37   1.49  ... -7.511  0.021   
63         1.37   1.49  ... -7.511  0.021   
64         1.37   1.49  ... -7.511  0.021   
65         1.37   1.49  ... -7.511  0.021   
66         1.37   1.49  ... -7.511  0.021   

                                            reference  couple_lower  \
0   N D        Kurucz Fe I 2014Fe               1 ...            LS   
1   E 0.05     Wisconsin REE ex(153)Eu+         3 ...            LS   
2   E 0.05     Wisconsin REE ex(151)Eu+         3 ...            LS   
3   E 0.05     Wisconsin REE ex(153)Eu+         3 ...            LS   
4   E 0.05     Wisconsin REE ex(151)Eu+         3 ...            LS   
..                                                ...           ...   
62  NC         Kurucz MnI 2007 Mn               8 ...            LS   
63  NC         Kurucz MnI 2007 Mn               8 ...            LS   
64  NC         Kurucz MnI 2007 Mn               8 ...            LS   
65  NC         Kurucz MnI 2007 Mn               8 ...            LS   
66  NC         Kurucz MnI 2007 Mn               8 ...            LS   

          term_lower couple_upper                       term_upper error  \
0            3d8 c3F           LS       3d6.(3F2).4s.4p.(3P*) v3D*  0.50   
1   4f7.(8S).5d a9D*           JJ  4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00   
2   4f7.(8S).5d a9D*           JJ  4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00   
3   4f7.(8S).5d a9D*           JJ  4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00   
4   4f7.(8S).5d a9D*           JJ  4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00   
..               ...          ...                              ...   ...   
62       3d5.4s2 b4D           LS                 3d6.(5D).4p z4D*  0.25   
63       3d5.4s2 b4D           LS                 3d6.(5D).4p z4D*  0.25   
64       3d5.4s2 b4D           LS                 3d6.(5D).4p z4D*  0.25   
65       3d5.4s2 b4D           LS                 3d6.(5D).4p z4D*  0.25   
66       3d5.4s2 b4D           LS                 3d6.(5D).4p z4D*  0.25   

   atom_number  ionization  
0          1.0         1.0  
1          1.0         2.0  
2          1.0         2.0  
3          1.0         2.0  
4          1.0         2.0  
..         ...         ...  
62         1.0         1.0  
63         1.0         1.0  
64         1.0         1.0  
65         1.0         1.0  
66         1.0         1.0  

[67 rows x 22 columns]
```

Now we calculate the central depth and line range, in Teff, logg and monh of 5000, 4.0, 0:

```py
sme = SME_Structure()
sme.teff, sme.logg, sme.monh = 5000, 4.0, 0
sme.linelist = linelist

synthesizer = Synthesizer()
sme = synthesizer.update_cdr(sme)
```

Then `sme.linelist` becomes:

```py
   species     wlcent  gflog     excit  j_lo     e_upp  j_up  lande_lower  \
0     Fe 1  6436.4055 -2.460  4.186364   2.0  6.112128   1.0         0.68   
1     Eu 2  6437.6106 -1.242  1.319612   5.0  3.245008   5.0         1.73   
2     Eu 2  6437.6121 -1.280  1.319636   5.0  3.245027   5.0         1.73   
3     Eu 2  6437.6135 -2.473  1.319612   5.0  3.245008   5.0         1.73   
4     Eu 2  6437.6194 -2.511  1.319636   5.0  3.245027   5.0         1.73   
..     ...        ...    ...       ...   ...       ...   ...          ...   
62    Mn 1  6443.4689 -2.818  3.772307   1.5  5.695960   2.5         1.21   
63    Mn 1  6443.4726 -3.538  3.772307   1.5  5.695960   2.5         1.21   
64    Mn 1  6443.4821 -2.275  3.772307   1.5  5.695960   2.5         1.21   
65    Mn 1  6443.4887 -2.964  3.772307   1.5  5.695960   2.5         1.21   
66    Mn 1  6443.4938 -3.918  3.772307   1.5  5.695960   2.5         1.21   

    lande_upper  lande  ...  couple_lower        term_lower  couple_upper  \
0          0.56   0.73  ...            LS           3d8 c3F            LS   
1          1.79   1.76  ...            LS  4f7.(8S).5d a9D*            JJ   
2          1.79   1.76  ...            LS  4f7.(8S).5d a9D*            JJ   
3          1.79   1.76  ...            LS  4f7.(8S).5d a9D*            JJ   
4          1.79   1.76  ...            LS  4f7.(8S).5d a9D*            JJ   
..          ...    ...  ...           ...               ...           ...   
62         1.37   1.49  ...            LS       3d5.4s2 b4D            LS   
63         1.37   1.49  ...            LS       3d5.4s2 b4D            LS   
64         1.37   1.49  ...            LS       3d5.4s2 b4D            LS   
65         1.37   1.49  ...            LS       3d5.4s2 b4D            LS   
66         1.37   1.49  ...            LS       3d5.4s2 b4D            LS   

                         term_upper error atom_number ionization  \
0        3d6.(3F2).4s.4p.(3P*) v3D*  0.50         1.0        1.0   
1   4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00         1.0        2.0   
2   4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00         1.0        2.0   
3   4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00         1.0        2.0   
4   4f7.(8S<7/2>).6p<3/2> (7/2,3/2)  1.00         1.0        2.0   
..                              ...   ...         ...        ...   
62                 3d6.(5D).4p z4D*  0.25         1.0        1.0   
63                 3d6.(5D).4p z4D*  0.25         1.0        1.0   
64                 3d6.(5D).4p z4D*  0.25         1.0        1.0   
65                 3d6.(5D).4p z4D*  0.25         1.0        1.0   
66                 3d6.(5D).4p z4D*  0.25         1.0        1.0   

   central_depth line_range_s  line_range_e  
0       0.796796    6435.9555     6436.8555  
1       0.250884    6437.3106     6437.9106  
2       0.234380    6437.3121     6437.9121  
3       0.018652    6437.3135     6437.9135  
4       0.017111    6437.3194     6437.9194  
..           ...          ...           ...  
62      0.149999    6443.1689     6443.7689  
63      0.031953    6443.1726     6443.7726  
64      0.388658    6443.1821     6443.7821  
65      0.111368    6443.1887     6443.7887  
66      0.013534    6443.1938     6443.7938  

[67 rows x 25 columns]
```

The last three columns are the central depth, start and end wavelength of line range.
The function calculated these parameters for each line in a batch of 2000 lines (`sme.cdr_N_line_chunk`), in parallel (`sme.cdr_parallel`) of 10 (`sme.cdr_n_jobs`) jobs. 
These vairalbes can be modified to change the behaviour of the function.

The variable `sme.linelist.cdr_paras` will also be updated to an array of [Teff, logg, monh, vmic] used for running `update_cdr`.

## How to synthesize a long spectra

The synthesis of long spectra, i.e., a spectra covering a wide wavelength range (>100AA) or a spectra divided into many chunks but still cover a wide wavelength range, would be a bit tricky.

The deafult synthesis mode of PySME takes all the lines in the `linelist` into account to synthesize the spectra, thus the synthesis takes a long time.
To make it even worse, the synthesis is repeated for each segment, thus if we have 10 segment covering from 3000-9000AA, all the lines (including those in ~9000AA) will be used for the first segment (3000-3600AA), and the total time cost is 10 times longer than just using one segment to cover the whole range.

This situation can be improved by using the central depth and line range parameters to select the relevant lines for each segment, by simply adding `linelist_mode='dynamic'` in `synthesize_spectrum`:

```py
from pysme.synthesize import synthesize_spectrum

sme = SME_Structure()
sme.teff, sme.logg, sme.monh, sme.vmic, sme.vmac, sme.vsini = 5777, 4.0, 0.07, 1.86, 0, 9.61
sme.wave = np.arange(3000, 6000, 0.02).reshape(-1, 100)
sme.linelist = line_list
sme.nlte.set_nlte('Na')
sme.cdr_depth_thres = 0.1
sme = synthesize_spectrum(sme, linelist_mode='dynamic')
```

`linelist_mode='auto'` is still accepted as a deprecated compatibility alias
and maps to `linelist_mode='dynamic'`.

The detailed explanaiton of what `linelist_mode='dynamic'` triggers is as follow:

1. If the `sme.linelist.cdr_paras` is None (means `update_cdr` hasn't run for it), or either the Teff, logg, monh and vmic (optional) in `sme.linelist.cdr_paras` differ from the input parameters by 250K, 0.5, 0.5 or 1 (you can change these parameter in `sme.linelist.cdr_paras_thres` as a dictionary), then run `update_cdr` under current input stellar parameters.
2. Find the beginning and ending wavelength of each segment, wbeg and wend. 
3. For each segment, only input the lines with:
    - `central_depth` > sme.cdr_depth_thres (default 0), and 
    - `line_range_s` <= wend + del_wav + line_margin, and 
    - `line_range_e` >= wbeg - del_wav - line_margin
4. Finish the synthesis.

Here del_wav is defined as sqrt(sme.vmic^2 + sme.vmac^2 + sme.vsini^2)/c*lambda + lambda/sme.ipres, and line_margin is set to 2AA (just for adding some margin).
NLTE calculation is also availalbe in this mode, by setting `sme.nlte.set_nlte`.

Note that the input line list will not change - only the lines being input to each segment synthesis will differ.  

Removing the `linelist_mode` variable or changing it to `all` falls back to the default way to perform synthesis, described in the beginning of this section.

## Unified line-selection controls (CDR and ALMAX)

Current PySME exposes a unified interface for external line selection:

- `sme.line_select_method`: `internal | cdr | almax`
- `sme.line_select_parallel`: enable/disable parallel metadata update
- `sme.line_select_n_jobs`: worker count for parallel updates
- `sme.line_select_chunk_size`: chunk size for metadata updates
- `sme.line_select_recompute`: `if_stale | always | never`

ALMAX mode supports two strong-line rules:

- `sme.line_select_almax_use_bins = False`:
  - simple cut with `almax_ratio >= sme.line_select_almax_threshold`
- `sme.line_select_almax_use_bins = True`:
  - bin-wise cumulative cut with `flag_strong_lines_by_bins(...)`

Both ALMAX rules use the same threshold parameter:

- `sme.line_select_almax_threshold`

If `sme.line_select_almax_threshold` is `None`, it falls back to `sme.accrt`
to preserve legacy behavior.

Example:

```py
sme.line_select_method = "almax"
sme.line_select_almax_threshold = sme.accrt
sme.line_select_almax_use_bins = True
sme.line_select_almax_bin_width = 0.2

sme = synthesize_spectrum(sme, linelist_mode="dynamic")
```

## How to fit parameter for a long spectra

The vairalbe for such process is as same as the one used in the last section:

```py
from pysme.solve import solve
sme = SME_Structure()
sme.teff, sme.logg, sme.monh, sme.vmic, sme.vmac, sme.vsini = 5777, 4.0, 0.07, 1.86, 0, 9.61
sme.wave = obs_wave
sme.spec = obs_flux
sme.uncs = obs_flux_err
sme.linelist = line_list
sme.cdr_depth_thres = 0.1
sme = solve(sme, linelist_mode='dynamic')
```

It only triggers the `linelist_mode` variable in `synthesize_spectrum`.
