Zemax(Optics Studio) を Python で制御するメモ

## 背景

Zemax(Optics Studio)を Python で制御したい...

C++ API もありますが, そこまで C++ を駆使するほどでもないので Python でぱぱっとやりたい

## 仕組み

大元は C# あたりで記載されているようで, 各言語は .NET binding で提供されているようです(すくなくとも, Python はそうなっている)

スタンドアローン(dll 直接叩いて処理. Optics Studio 起動しなくとも動く),
Optics Studio の画面に結果を出すなどのインタラクティブモードの２つのモードがあります.

https://support.zemax.com/hc/en-us/articles/1500005489681-Connect-to-OpticStudio-from-Python-with-the-ZOS-API-Interactive-Extension

ここではスタンドアローンでやります.
Python であれば Matplotlib などでプロットしたりもできますしね.

また, インタラクティブモードでもある程度動きますが, Optics Studio の GUI との連携はとれません...
Optics Studio 側で Surface の値など編集して, それを interactive にPython でなにか計算させるなどはできない. サポートからの反応も「なんか interactive で連携するには内部構造からしていろいろ難しいぽよ...」とありました... ので, もしシームレスに対応されるとしても, たぶん数年は無理っぽそうです.


## Python 3.8?

現時点(2022/05/10)の最新 22.1.2 では, Python 3.8 未対応(PythonNET が未対応)とありますが,
とりあえず Python3.8 でも pip で入れればいけます.

```
$ python -m pip install pythonet==2.5.2
```

最新 3.x バージョンではいくつかおかしい(メソッド呼ぶと落ちたりする)のを確認しました.
Python 3.8 + 2.5.2 で動作確認いたしました.
(Pyhthon 3.7 だと pythonnet 最新で動くかもしれません)

## JupyerLab desktop

いろいろ API 試すのであれば, JupyterLab desktop でまずはいろいろやって見るのがよいでしょう!

Jupyterlab の場合は

```
%pip install pythonnet
```

なりして `pythonnet` 入れておきます.


jupyterlab-lsp を入れておくと zos モジュールのメソッドが補完できて便利です!

https://yiskw713.hatenablog.com/entry/jupyterlab-setup

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">人生初 JupyterLab + Python lsp(便利すぎた...)✌️ <a href="https://t.co/nej9TWkCEb">pic.twitter.com/nej9TWkCEb</a></p>&mdash; Syoyo Fujita 🌸 レイトラ ® 🐯 4 周年 🎉 (@syoyo) <a href="https://twitter.com/syoyo/status/1524660369721262080?ref_src=twsrc%5Etfw">May 12, 2022</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script> 


### 情報

https://github.com/xzos/PyZDDE

(ただし DDE ベースのやりかたはもうサポートされていないので, なんとなく仕組みチェックに利用するとよいでしょう)

https://support.zemax.com/hc/en-us/articles/1500005490181-Understanding-the-main-Interfaces-of-ZOS-API-using-MATLAB


### 構成

```
zos.TheApplication
zos.TheSystem
```

基本 `C#` へのバインディングなので, print とかしても内容うまく表示してくれないので注意しましょう...

### 無限 object thickness

float('inf')でいけるはずです.

OpticsStudio では `1.0e+10` 以上が無限として扱われますので, `1.0e+10` を指定でもよいかもしれません.

### ファイルを読み込む

`TheApplication.PrimarySystem.LoadFile`

ZAR(Archive 形式)を読み込むのもできます. API 探してみてくださいネ.

### (quick)フォーカシング

https://support.zemax.com/hc/en-us/articles/1500005489681-Connect-to-OpticStudio-from-Python-with-the-ZOS-API-Interactive-Extension

`TheSystem.Tools.OpenQuickFocus()`


### first-order (optics) data

日本語では近軸近似量?

ZOS API では Operand 経由で一個一個取得します. ちょっとめんどいですね.

たとえば EXPP(Exit pupil position)

```
TheMFE = TheSystem.MFE

OpEXPP = ZOSAPI.Editors.MFE.MeritOperandType.EXPP

TheMFE.GetOperandValue(OpEXPP, 0, 0, 0, 0, 0, 0, 0, 0)
```

### Field(像高)設定

`TheSystem.SystemData` で行います.
e.g.

```
TheSystemData = zos.TheSystem.SystemData
Field_1 = TheSystemData.Fields.GetField(1)
NewField_2 = TheSystemData.Fields.AddField(0, 5.0, 1.0)
```

最初の Field は必ず存在するので, 設定する場合は `AddField` ではなく,
`GetField(1).X = 0.0` のようにします.

`AddField` では `(x, y, weight)` だけしか指定できないので,
VDX などは `AddField` したのち, `GetField(2).VDX = 0.01` のようにして設定します.


`GetFieldType()` で field type を取得, `SetFieldType` で field type 設定できます.
(カッコは OpTaliX に対応するもの)

```
0 : Angle(XAN, YAN)
1 : Object height(XOB, YOB)
2 : Paraxial image height(XIM, YIM)
3 : Real image height(XRI, YRI)
4 : Theodolite angle (天文用)
``` 


Object height の場合, レンズ(Surface editor)で無限遠を設定(Object plane を Infinity)できませんので注意です!(エラーが出る)

https://support.zemax.com/hc/en-us/articles/1500005576302-How-to-design-a-singlet-lens-Part-1-Setup

API で, Vignetting の計算とクリアもできます.

#### Field データの切り替え 

システム上では複数の Field テーブルセットを持つのはできないようです.
また, フィールドファイル `.fld` を直接読む API は無いっぽい?

`DeleteAllFields()` ですべての像高を削除し, 再度登録になるでしょう.
`DeleteAllFields()` ではリターンでいくつ消したか帰ってきます
(Field 1 個はかならず生成(残る)されるので, すくなくとも 1 以上が帰ってくる)


### 波長の設定

![Screenshot from 2022-07-06 17-17-47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/36275/1cd367d7-a3a2-8841-ae85-d043775d15ed.png)


Field と同じようにして,
`TheSystem.Wavelengths` で指定.

波長データは 24 個固定で, `DeleteAllWavelengths` みたいなメソッドはありません.
設定し直ししたいときは,

`RemoveWavelength(1)` のようにして該当の wavelength 情報を消してから, AddWavelength(これはしかし番号指定はできず, 開いている番号に追加されるっぽい)

API ではどれが active(GUI でチェックが入っているか)かは `.IsActive` でわかります.
`MakeActive` みたいなのは無いようです.

`AddWavelength(wavelen, weight)` で波長のデータを追加できます.

ちなみに単位は `[um]` 固定のようです..

波長にはいくつか Preset が用意されています.

`ZOSAPI.SystemData.WavelengthPreset.FdC_Visible` など.

Preset を使う場合,

`SelectWavelengthPreset(preset)` を使います.

こちらも, 波長ファイル `.wav` を直接読む API はなさそうです...
 
#### Primary index の設定

`GetWavelength(n).MakePrimary()` で設定できます.

`.IsPrimary` で Primary かどうかわかります.

`NumberOfWavelengths` で(アクティブな)波長の数がわかりますが,
どの wavelength が primary かはアクティブな波長をスキャンして見つけるしかななさそうです.



### FFT PSF を計算する

https://community.zemax.com/zos-api-12/how-can-i-set-imagedelta-of-fftpsf-in-zosapi-matlab-1155

https://community.zemax.com/zos-api-12/settings-syntax-for-python-535

DataGrid に結果があります.

データ自体は .NET でのデータなので, numpy 配列にするには `zos.reshape` あたりで変換します.

```py
TheAnalyses = TheSystem.Analyses
psf = TheAnalyses.New_FftPsf()

psf_settings = newWin.GetSettings()
psf_settings.SampleSize = ZOSAPI.Analysis.SampleSizes.S_256x256
...

# TODO: 
# compute
psf.ApplyAndWaitForCompletion()

ret = psf.GetResults()
print(ret.NumberOfDataGrids)

grid = ret.GetDataGrid(0)

raw = grid.ValueData().get_Data()
print(raw.GetLength(0), raw.GetLength(1))

m = zos.reshape(raw, raw.GetLength(0), raw.GetLength(1)


import matplotlib.pyplot as plt
plt.imshow(m)
```

![Screenshot from 2022-05-12 21-08-19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/36275/f7f097f2-944c-f477-e0d9-699e4ae6f8ff.png)

Voila!

PSF の場合は(複数 Fields 一括処理や Grid PSF 計算でなければ?) Result の NumberOfDataGrids は 1 になります.


#### 設定

 ![Screenshot from 2022-06-27 18-09-56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/36275/689795e2-e077-aa6f-a24f-8c101e98e41a.png)

https://community.zemax.com/zos-api-12/how-can-i-set-imagedelta-of-fftpsf-in-zosapi-matlab-1155

* SampleSize(Sampling)
* OutputSize(Display)
* ImageDelta(double value)
  * 像面(image plane)でのサンプル点の間隔を設定. 0 = デフォルトを使う. 
  * ImageDelta * SampleSize が PSF の計算領域となる https://support.zemax.com/hc/ja/articles/1500005576562-%E7%82%B9%E5%83%8F%E5%BC%B7%E5%BA%A6%E5%88%86%E5%B8%83%E9%96%A2%E6%95%B0%E3%81%A8%E3%81%AF
  * negative 値を取ることもできる(Zemax manual 参照)
  * Optalix での patch size と同等かと思われる. 計算式などは OpTaliX reference manual http://www.optenso.com/download/download.html  14.2.2.1 も参考になるでしょう.
* Field
  * Field 番号(像高の番号)を指定します. `Field` は `IAS_Field` オブジェクトのため, `settings.Field.SetFieldNumber(1)` のようにしてフィールド番号を指定します.
  * `UseAllFields` もありますが, これを指定すると全像高で一括 PSF 計算してくれるのかしらん?(未検証)
* Wavelength
  * こちらも Field と同じように `IAS_Wavelength` オブジェクトです.
  * `settings.Wavelength.SetWavelengthNumber(1)` や, 全波長を使う(デフォルト)場合は `settings.Wavelength.UseAllWavelengths()` とします. 

FFT PSF の場合, OutputSize(Display) は SampleSize の最大 2 倍までにしないとうまく計算されないので注意です!
(e.g. SampleSize 64x64 なら, OutputSize は 128x128 まで)



#### Huygens PSF

![Screenshot from 2022-06-27 18-18-35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/36275/abdd89a9-356b-523e-fad3-40617e7aed65.png)

`New_HuygensPsf` を使います.

* PupilSampleSize(Pupil Sampling)
* ImageSampleSize(Image Sampling)
* ImageDelta
  * 0.0 の場合は自動で計算される. この自動で計算のときの式は FFT PSF のとは異るので注意 



#### 計算のステートや中断

Huygens PSF などは, 解像度高いと計算時間かかってしまします.
処理の中断なども API で制御したい場合は,

* `Apply()` で計算を開始する
* `IsRunning()` で処理中かチェック
* `Close()` で処理が終わっていなければ中断

とします.

https://community.zemax.com/zos-api-12/zos-api-unexpectedly-fails-to-open-the-wavefront-analysis-window-2130

`Apply()` および `ApplyAndWaitForCompletion()` では, 現在計算中のものは terminate されます.

#### Strehl ratio

`STRH` Operand で取得します.

sampling level(1 = 32x32, 2 = 64x64, ...), wavelength index(0 だと全波長), field を指定します
(polarization, `All conf`(?)も必要であれば)

![Screenshot from 2022-07-07 18-57-31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/36275/994430eb-60c8-ff91-022b-3d1a7ef38d41.png)

```
TheMFE = TheSystem.MFE

OpSTRH = ZOSAPI.Editors.MFE.MeritOperandType.STRH

samp = ...
wave = ...
field = ...

TheMFE.GetOperandValue(OpSTRH, samp, wave, field, 0, 0, 0, 0, 0)
```

### MTF を計算する

MTF はサンプルスクリプトに例があるのでそちらをみるのがいいでしょう.
DataSeries に結果があります.

ちなみに Through Focus では image plane(sensor image plane) を 0 基準に計算されます
(paraxial image plane ではない)

古い Zemax では Defocus 基準(paraxial image plane 基準)だったような? :thinking:

### Analysis

PSF, MTF とか, このあたりはまとめて Analysis オブジェクトとして管理されています.

* Title でタイトル取得
* GetAnalysisName で Analysis の名称

が取得できます. Title と AnalysisName はだいたい似たような文字列が得られます(たぶん).

#### Text ファイルの取得

```
ret = analysis.GetResults()

ret.GetTextFile(filename)
```

で `.txt` ファイルが取得できます.

ファイルパスは相対パスだと出力されないようです(return status は true であるが... Zemax のデフォルト? フォルダに出力されているのかしら?).

`GetTextFile(os.path.join(os.getcwd(), "bora_psf.txt"))` のようにして絶対パス指定してあげましょう.

また, デフォルトだと UTF-16LE(Windows Unicode)のエンコードになります.
(Zemax 側の設定で Ascii(CP932)にすることもできるが, なんかあんまりもうサポートしていないようで文字化けしたりする)

### Optimization

Zemax での最適化はオペランドで最適化条件(制約条件)を記述していって最適化するため, プログラマブルな感じです.
(たとえば OpTaliX のようにぺろっと MTF で autofocus させるとかはできない)
一応 Quickfocus 機能ありますが, これは最適化条件を設定するためのプリセットみたいな感じです.

Autofocus するには, レンズ最背面 Thickness を変数にし, この Field の MTF 値(MTF も FFT MTF(MTFA/S/T) や Huygens MTF(MTHA/S/T) などオペランドがいくつかある)が最小になるのを探す.. みたいな感じになります.
 
#### Optimization wizard

`TheSystem.TheMFE.SEQOptimizationWizard2` にあります.
最近の Optics Studio での GUI ででる Optimization Wizard はこの `2` になります.

その他最適化最適化は sample にありますのでそれを参考にするとよいでしょう.

また, OptimizationWizard では各種オペランドを一括で生成するだけの感じです.


#### 境界条件

オペランド設定してがんばります..

#### Thickness 境界の設定

たとえば autofocus などで最背面の thickness を調整したいが, それに可動範囲を設定したいとき.

http://www.lensya.co.jp/010/wforum.cgi?mode=allread&no=2129&page=0

CTGT
CTLT 

でいけます.

最適化 wizard の設定で生成される MXCA あたりでは, air か glass かも考慮してより細かく制御できるようです(未検証)

#### Configuration

Optimization Wizard など, Configuration を指定することができます.
Configuration はいくつかの設定をまとめて, 切り替える感じの機能です. 
最適化のときや PSF 計算のときに設定切り替えが楽になります.

#### コントラスト最適化

Autofucus(MTF 最適化)の場合はコントラスト最適化がよいようです.

https://support.zemax.com/hc/ja/articles/1500005488961-%E3%82%B3%E3%83%B3%E3%83%88%E3%83%A9%E3%82%B9%E3%83%88%E6%9C%80%E9%81%A9%E5%8C%96%E3%82%92%E4%BD%BF%E7%94%A8%E3%81%97%E3%81%9F-MTF-%E6%80%A7%E8%83%BD%E3%81%AE%E6%9C%80%E9%81%A9%E5%8C%96

ZOS Python では以下のような感じにしたらあとはよろしく operand ぺろっと生成してくれます.

```py
    TheMFE = TheSystem.MFE
    OptWizard = TheMFE.SEQOptimizationWizard2 # Must be 2
    OptWizard.Criterion = self.ZOSAPI.Wizards.CriterionTypes.Contrast
    OptWizard.SpatialFrequency = 30 # lp/mm
```

Exit pupil 面に Rings でサンプル点(hexapolar かしらん? https://qiita.com/syoyo/items/4fffee1d460d14f9336c )を生成し, レイトレして波面収差(OPD)を最小にする感じですかね.
その値を最適化するようなオペランドが生成されます.


### レンズデータへのアクセス

LTE(Lens Data Editor)を使います.

`GetSurfaceAt(i)` で i 番目の surface を取得できます.
 
```py
TheLDE = TheSystem.LDE 

s0 = TheLDE.GetSurfaceAt(0)
```

* TypeName : Surface type name. `Standard` など
  * Type ですと id が得られる.
* Radius: Radius
* Thickness: Thickness
* Material : ガラス(文字列で指定/取得)
* Conic: Conic 
* Comment : コメント
* IsObject: Object plane かどうか
* IsStop: Stop(Aperture) かどうか
* IsImage: Image plane かどうか

#### 非球面係数

Surface type により, メソッドが異なります.

Even Aspherical の場合

```
asph_s = TheLDE.GetSurfaceAt(n)
asph_s.SurfaceData.GetNthEvenOrderTerm(4)
```

で取得できます! N 次は 2, 4, 6, ...

### Glass catalog へのアクセス

`TheSystem.SystemData.MaterialCatalogs`


`GetMaterialCatalogs()` でリストを取得できます.

```
mats = TheSystem.SystemData.MaterialCatalogs.GetMaterialCatalogs()
for m in mats:
  print(m)
```

`GetMaterialInCatalog(catalog)` で Catalog から Material のリスト取得できます.

ただ, これですと単に Material の文字列(e.g. `BK7`)を得られるだけです.


### Material 情報へのアクセス

残念ながら Material のデータ構造(屈折率や, Sellmeier 係数とか)には ZOS API からはアクセスできません :cry:
値がほしいときは Glass file を直読みするか, Macro かなんかオペランド指定して評価したら取得できるかもしれません.


カスタムの Material(Glass) を扱いたい場合, GUI でやるか, Glass file をいじってそれを与えるしかなさそうです.

#### Glass file format

AGF : Ascii 形式(BGF : Binary 形式(Zemax が AGF から自動で作る))
ZTG : Table(index, value のペアのリスト)

https://support.zemax.com/hc/en-us/articles/1500005486681-Zemax-file-extensions

その他 GRD とかもありますが, 基本は AGF, ZTG でしょうか.

#### GRIN material

T.B.W.

### レイトレーシング(光線追跡)

T.B.W.

#### Volume scattering

Plugin? を使うことで subsurface scattering の計算ができました :tada:

https://support.zemax.com/hc/ja/articles/1500005579202?input_string=get+optical+path+ray+data+in+zos-api%28sequential+mode%29



### Zemax macro を Python から呼びたい

ZOS API ではなくて, Zemax macros を叩きたいときもあります.

https://community.zemax.com/zos-api-12/run-a-macro-from-zosapi-769

や PyZDDE を参照してみましょう.

## Operand

よくつかいそうなの

* EFFL : Effective Focal Length (波長指定)
* ISFN : Image Space Fno (いわゆる F 値)
* MTFA : FFT MTF average((sagittal + tangential)/2)
* MTHA : Huygens MTF average((sagittal + tangential)/2)

## ファイルエンコーディング

Zemax では BOM 付き utf16-le で .TXT など書き出されます.
一応 Zemax 側の設定で Ascii(日本語環境だと CP932)で出す設定もありますが, 表示がバグったりするのでおすすめしません.

python の場合は open 時に `encoding="utf-16"` を指定してあげます.

zemax .txt を編集して再度 zemax に読ませる場合(e.g. 波長データなど), 再度 utf16-le でエンコードするか,

https://qiita.com/syoyo/items/8d0d054eabc127584ffe

iconv で変換してあげる必要があります.