From d6b636af63071b9000854edc92862df76dfdcd72 Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 17:02:38 -0400 Subject: [PATCH 1/7] feat: prepare tested plane-wave optics preview Package the bounded plane-wave tool with numerical benchmarks, buyer exports, license notices and draft-only release automation. Preserve explicit scientific and commercial publication gates. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .github/workflows/plane-wave-preview.yml | 117 +++++++ .gitignore | 11 + LICENSE | 202 ++++++++++++ LICENSE-CC-BY-4.0.txt | 396 +++++++++++++++++++++++ NOTICE | 18 ++ README.md | 118 ++++++- pyproject.toml | 39 +++ tests/test_tmm_oled9.py | 213 ++++++++++++ tmm_OLED9.py | 254 ++++++++++----- 9 files changed, 1287 insertions(+), 81 deletions(-) create mode 100644 .github/workflows/plane-wave-preview.yml create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 LICENSE-CC-BY-4.0.txt create mode 100644 NOTICE create mode 100644 pyproject.toml create mode 100644 tests/test_tmm_oled9.py diff --git a/.github/workflows/plane-wave-preview.yml b/.github/workflows/plane-wave-preview.yml new file mode 100644 index 0000000..53d92cd --- /dev/null +++ b/.github/workflows/plane-wave-preview.yml @@ -0,0 +1,117 @@ +name: Plane-wave preview + +on: + push: + branches: [main] + pull_request: + workflow_dispatch: + inputs: + draft_tag: + description: "Optional existing tag, e.g. v0.1.0a1, for a draft prerelease only" + required: false + type: string + +permissions: + contents: read + +jobs: + test: + runs-on: ${{ matrix.os }} + timeout-minutes: 15 + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, windows-latest] + python: ['3.11', '3.12', '3.13', '3.14'] + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python }} + - run: python -m pip install ".[dev]" + - run: python -m pytest -q --junitxml=test-results.xml + - uses: actions/upload-artifact@v4 + if: always() + with: + name: results-${{ matrix.os }}-${{ matrix.python }} + path: test-results.xml + if-no-files-found: error + + package: + needs: test + runs-on: ubuntu-latest + timeout-minutes: 15 + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - run: python -m pip install ".[dev]" + - run: python -m build + - run: python -m twine check dist/* + - name: Test installed wheel outside source tree + shell: bash + run: | + python -m venv "$RUNNER_TEMP/buyer" + "$RUNNER_TEMP/buyer/bin/python" -m pip install dist/*.whl + cd "$RUNNER_TEMP" + "$RUNNER_TEMP/buyer/bin/mongelli-plane-wave" --wavelength-count 3 --angle-count 3 --output preview.png + test -s preview.png && test -s preview.csv && test -s preview.json + - name: Audit installed runtime dependencies and produce CycloneDX inventory + run: | + SITE=$("$RUNNER_TEMP/buyer/bin/python" -c "import sysconfig; print(sysconfig.get_paths()['purelib'])") + python -m pip_audit --path "$SITE" --format cyclonedx-json --output dist/dependency-audit.cdx.json + - name: Record artifact checksums + run: sha256sum dist/*.whl dist/*.tar.gz dist/*.json > dist/SHA256SUMS.txt + - uses: actions/upload-artifact@v4 + with: + name: plane-wave-preview-distributions + path: dist/* + if-no-files-found: error + + draft: + if: github.event_name == 'workflow_dispatch' && inputs.draft_tag != '' + needs: package + runs-on: ubuntu-latest + timeout-minutes: 5 + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - name: Validate version and existing tag targets the tested commit + env: + DRAFT_TAG: ${{ inputs.draft_tag }} + run: | + python - <<'PY' + import os + import subprocess + import tomllib + from pathlib import Path + version = tomllib.loads(Path("pyproject.toml").read_text())["project"]["version"] + tag = os.environ["DRAFT_TAG"] + if tag != "v" + version: + raise SystemExit("Tag must match package version") + target = subprocess.check_output(["git", "rev-parse", "--verify", f"refs/tags/{tag}^{{commit}}"], text=True).strip() + head = subprocess.check_output(["git", "rev-parse", "HEAD"], text=True).strip() + if target != head: + raise SystemExit("Tag must target the tested workflow commit") + PY + - uses: actions/download-artifact@v4 + with: + name: plane-wave-preview-distributions + path: dist + - name: Create draft only + env: + GH_TOKEN: ${{ github.token }} + DRAFT_TAG: ${{ inputs.draft_tag }} + run: | + gh release create "$DRAFT_TAG" dist/* --verify-tag --draft --prerelease \ + --title "Plane-wave optics preview $DRAFT_TAG" \ + --notes "Illustrative plane-wave transmittance only. Not a validated OLED dipole-outcoupling solver. Not approved for paid publication; see README commercial approval gates." diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ad9efa1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +.venv/ +.pytest_cache/ +__pycache__/ +*.pyc +*.egg-info/ +build/ +dist/ +.serena/ +transmission.png +transmission.csv +transmission.json diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d645695 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/LICENSE-CC-BY-4.0.txt b/LICENSE-CC-BY-4.0.txt new file mode 100644 index 0000000..da6ab6c --- /dev/null +++ b/LICENSE-CC-BY-4.0.txt @@ -0,0 +1,396 @@ +Attribution 4.0 International + +======================================================================= + +Creative Commons Corporation ("Creative Commons") is not a law firm and +does not provide legal services or legal advice. Distribution of +Creative Commons public licenses does not create a lawyer-client or +other relationship. Creative Commons makes its licenses and related +information available on an "as-is" basis. Creative Commons gives no +warranties regarding its licenses, any material licensed under their +terms and conditions, or any related information. Creative Commons +disclaims all liability for damages resulting from their use to the +fullest extent possible. + +Using Creative Commons Public Licenses + +Creative Commons public licenses provide a standard set of terms and +conditions that creators and other rights holders may use to share +original works of authorship and other material subject to copyright +and certain other rights specified in the public license below. The +following considerations are for informational purposes only, are not +exhaustive, and do not form part of our licenses. + + Considerations for licensors: Our public licenses are + intended for use by those authorized to give the public + permission to use material in ways otherwise restricted by + copyright and certain other rights. Our licenses are + irrevocable. Licensors should read and understand the terms + and conditions of the license they choose before applying it. + Licensors should also secure all rights necessary before + applying our licenses so that the public can reuse the + material as expected. Licensors should clearly mark any + material not subject to the license. This includes other CC- + licensed material, or material used under an exception or + limitation to copyright. More considerations for licensors: + wiki.creativecommons.org/Considerations_for_licensors + + Considerations for the public: By using one of our public + licenses, a licensor grants the public permission to use the + licensed material under specified terms and conditions. If + the licensor's permission is not necessary for any reason--for + example, because of any applicable exception or limitation to + copyright--then that use is not regulated by the license. Our + licenses grant only permissions under copyright and certain + other rights that a licensor has authority to grant. Use of + the licensed material may still be restricted for other + reasons, including because others have copyright or other + rights in the material. A licensor may make special requests, + such as asking that all changes be marked or described. + Although not required by our licenses, you are encouraged to + respect those requests where reasonable. More considerations + for the public: + wiki.creativecommons.org/Considerations_for_licensees + +======================================================================= + +Creative Commons Attribution 4.0 International Public License + +By exercising the Licensed Rights (defined below), You accept and agree +to be bound by the terms and conditions of this Creative Commons +Attribution 4.0 International Public License ("Public License"). To the +extent this Public License may be interpreted as a contract, You are +granted the Licensed Rights in consideration of Your acceptance of +these terms and conditions, and the Licensor grants You such rights in +consideration of benefits the Licensor receives from making the +Licensed Material available under these terms and conditions. + + +Section 1 -- Definitions. + + a. Adapted Material means material subject to Copyright and Similar + Rights that is derived from or based upon the Licensed Material + and in which the Licensed Material is translated, altered, + arranged, transformed, or otherwise modified in a manner requiring + permission under the Copyright and Similar Rights held by the + Licensor. For purposes of this Public License, where the Licensed + Material is a musical work, performance, or sound recording, + Adapted Material is always produced where the Licensed Material is + synched in timed relation with a moving image. + + b. Adapter's License means the license You apply to Your Copyright + and Similar Rights in Your contributions to Adapted Material in + accordance with the terms and conditions of this Public License. + + c. Copyright and Similar Rights means copyright and/or similar rights + closely related to copyright including, without limitation, + performance, broadcast, sound recording, and Sui Generis Database + Rights, without regard to how the rights are labeled or + categorized. For purposes of this Public License, the rights + specified in Section 2(b)(1)-(2) are not Copyright and Similar + Rights. + + d. Effective Technological Measures means those measures that, in the + absence of proper authority, may not be circumvented under laws + fulfilling obligations under Article 11 of the WIPO Copyright + Treaty adopted on December 20, 1996, and/or similar international + agreements. + + e. Exceptions and Limitations means fair use, fair dealing, and/or + any other exception or limitation to Copyright and Similar Rights + that applies to Your use of the Licensed Material. + + f. Licensed Material means the artistic or literary work, database, + or other material to which the Licensor applied this Public + License. + + g. Licensed Rights means the rights granted to You subject to the + terms and conditions of this Public License, which are limited to + all Copyright and Similar Rights that apply to Your use of the + Licensed Material and that the Licensor has authority to license. + + h. Licensor means the individual(s) or entity(ies) granting rights + under this Public License. + + i. Share means to provide material to the public by any means or + process that requires permission under the Licensed Rights, such + as reproduction, public display, public performance, distribution, + dissemination, communication, or importation, and to make material + available to the public including in ways that members of the + public may access the material from a place and at a time + individually chosen by them. + + j. Sui Generis Database Rights means rights other than copyright + resulting from Directive 96/9/EC of the European Parliament and of + the Council of 11 March 1996 on the legal protection of databases, + as amended and/or succeeded, as well as other essentially + equivalent rights anywhere in the world. + + k. You means the individual or entity exercising the Licensed Rights + under this Public License. Your has a corresponding meaning. + + +Section 2 -- Scope. + + a. License grant. + + 1. Subject to the terms and conditions of this Public License, + the Licensor hereby grants You a worldwide, royalty-free, + non-sublicensable, non-exclusive, irrevocable license to + exercise the Licensed Rights in the Licensed Material to: + + a. reproduce and Share the Licensed Material, in whole or + in part; and + + b. produce, reproduce, and Share Adapted Material. + + 2. Exceptions and Limitations. For the avoidance of doubt, where + Exceptions and Limitations apply to Your use, this Public + License does not apply, and You do not need to comply with + its terms and conditions. + + 3. Term. The term of this Public License is specified in Section + 6(a). + + 4. Media and formats; technical modifications allowed. The + Licensor authorizes You to exercise the Licensed Rights in + all media and formats whether now known or hereafter created, + and to make technical modifications necessary to do so. The + Licensor waives and/or agrees not to assert any right or + authority to forbid You from making technical modifications + necessary to exercise the Licensed Rights, including + technical modifications necessary to circumvent Effective + Technological Measures. For purposes of this Public License, + simply making modifications authorized by this Section 2(a) + (4) never produces Adapted Material. + + 5. Downstream recipients. + + a. Offer from the Licensor -- Licensed Material. Every + recipient of the Licensed Material automatically + receives an offer from the Licensor to exercise the + Licensed Rights under the terms and conditions of this + Public License. + + b. No downstream restrictions. You may not offer or impose + any additional or different terms or conditions on, or + apply any Effective Technological Measures to, the + Licensed Material if doing so restricts exercise of the + Licensed Rights by any recipient of the Licensed + Material. + + 6. No endorsement. Nothing in this Public License constitutes or + may be construed as permission to assert or imply that You + are, or that Your use of the Licensed Material is, connected + with, or sponsored, endorsed, or granted official status by, + the Licensor or others designated to receive attribution as + provided in Section 3(a)(1)(A)(i). + + b. Other rights. + + 1. Moral rights, such as the right of integrity, are not + licensed under this Public License, nor are publicity, + privacy, and/or other similar personality rights; however, to + the extent possible, the Licensor waives and/or agrees not to + assert any such rights held by the Licensor to the limited + extent necessary to allow You to exercise the Licensed + Rights, but not otherwise. + + 2. Patent and trademark rights are not licensed under this + Public License. + + 3. To the extent possible, the Licensor waives any right to + collect royalties from You for the exercise of the Licensed + Rights, whether directly or through a collecting society + under any voluntary or waivable statutory or compulsory + licensing scheme. In all other cases the Licensor expressly + reserves any right to collect such royalties. + + +Section 3 -- License Conditions. + +Your exercise of the Licensed Rights is expressly made subject to the +following conditions. + + a. Attribution. + + 1. If You Share the Licensed Material (including in modified + form), You must: + + a. retain the following if it is supplied by the Licensor + with the Licensed Material: + + i. identification of the creator(s) of the Licensed + Material and any others designated to receive + attribution, in any reasonable manner requested by + the Licensor (including by pseudonym if + designated); + + ii. a copyright notice; + + iii. a notice that refers to this Public License; + + iv. a notice that refers to the disclaimer of + warranties; + + v. a URI or hyperlink to the Licensed Material to the + extent reasonably practicable; + + b. indicate if You modified the Licensed Material and + retain an indication of any previous modifications; and + + c. indicate the Licensed Material is licensed under this + Public License, and include the text of, or the URI or + hyperlink to, this Public License. + + 2. You may satisfy the conditions in Section 3(a)(1) in any + reasonable manner based on the medium, means, and context in + which You Share the Licensed Material. For example, it may be + reasonable to satisfy the conditions by providing a URI or + hyperlink to a resource that includes the required + information. + + 3. If requested by the Licensor, You must remove any of the + information required by Section 3(a)(1)(A) to the extent + reasonably practicable. + + 4. If You Share Adapted Material You produce, the Adapter's + License You apply must not prevent recipients of the Adapted + Material from complying with this Public License. + + +Section 4 -- Sui Generis Database Rights. + +Where the Licensed Rights include Sui Generis Database Rights that +apply to Your use of the Licensed Material: + + a. for the avoidance of doubt, Section 2(a)(1) grants You the right + to extract, reuse, reproduce, and Share all or a substantial + portion of the contents of the database; + + b. if You include all or a substantial portion of the database + contents in a database in which You have Sui Generis Database + Rights, then the database in which You have Sui Generis Database + Rights (but not its individual contents) is Adapted Material; and + + c. You must comply with the conditions in Section 3(a) if You Share + all or a substantial portion of the contents of the database. + +For the avoidance of doubt, this Section 4 supplements and does not +replace Your obligations under this Public License where the Licensed +Rights include other Copyright and Similar Rights. + + +Section 5 -- Disclaimer of Warranties and Limitation of Liability. + + a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE + EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS + AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF + ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS, + IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION, + WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR + PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS, + ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT + KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT + ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU. + + b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE + TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION, + NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT, + INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES, + COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR + USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN + ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR + DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR + IN PART, THIS LIMITATION MAY NOT APPLY TO YOU. + + c. The disclaimer of warranties and limitation of liability provided + above shall be interpreted in a manner that, to the extent + possible, most closely approximates an absolute disclaimer and + waiver of all liability. + + +Section 6 -- Term and Termination. + + a. This Public License applies for the term of the Copyright and + Similar Rights licensed here. However, if You fail to comply with + this Public License, then Your rights under this Public License + terminate automatically. + + b. Where Your right to use the Licensed Material has terminated under + Section 6(a), it reinstates: + + 1. automatically as of the date the violation is cured, provided + it is cured within 30 days of Your discovery of the + violation; or + + 2. upon express reinstatement by the Licensor. + + For the avoidance of doubt, this Section 6(b) does not affect any + right the Licensor may have to seek remedies for Your violations + of this Public License. + + c. For the avoidance of doubt, the Licensor may also offer the + Licensed Material under separate terms or conditions or stop + distributing the Licensed Material at any time; however, doing so + will not terminate this Public License. + + d. Sections 1, 5, 6, 7, and 8 survive termination of this Public + License. + + +Section 7 -- Other Terms and Conditions. + + a. The Licensor shall not be bound by any additional or different + terms or conditions communicated by You unless expressly agreed. + + b. Any arrangements, understandings, or agreements regarding the + Licensed Material not stated herein are separate from and + independent of the terms and conditions of this Public License. + + +Section 8 -- Interpretation. + + a. For the avoidance of doubt, this Public License does not, and + shall not be interpreted to, reduce, limit, restrict, or impose + conditions on any use of the Licensed Material that could lawfully + be made without permission under this Public License. + + b. To the extent possible, if any provision of this Public License is + deemed unenforceable, it shall be automatically reformed to the + minimum extent necessary to make it enforceable. If the provision + cannot be reformed, it shall be severed from this Public License + without affecting the enforceability of the remaining terms and + conditions. + + c. No term or condition of this Public License will be waived and no + failure to comply consented to unless expressly agreed to by the + Licensor. + + d. Nothing in this Public License constitutes or may be interpreted + as a limitation upon, or waiver of, any privileges and immunities + that apply to the Licensor or You, including from the legal + processes of any jurisdiction or authority. + + +======================================================================= + +Creative Commons is not a party to its public +licenses. Notwithstanding, Creative Commons may elect to apply one of +its public licenses to material it publishes and in those instances +will be considered the “Licensor.” The text of the Creative Commons +public licenses is dedicated to the public domain under the CC0 Public +Domain Dedication. Except for the limited purpose of indicating that +material is shared under a Creative Commons public license or as +otherwise permitted by the Creative Commons policies published at +creativecommons.org/policies, Creative Commons does not authorize the +use of the trademark "Creative Commons" or any other trademark or logo +of Creative Commons without its prior written consent including, +without limitation, in connection with any unauthorized modifications +to any of its public licenses or any other arrangements, +understandings, or agreements concerning use of licensed material. For +the avoidance of doubt, this paragraph does not form part of the +public licenses. + +Creative Commons may be contacted at creativecommons.org. + diff --git a/NOTICE b/NOTICE new file mode 100644 index 0000000..0786d76 --- /dev/null +++ b/NOTICE @@ -0,0 +1,18 @@ +Mongelli plane-wave preview +Copyright (c) 2025 Guy Francis Mongelli + +Project code is Apache-2.0 under the existing repository notice. +Documentation, prose and figures are CC BY 4.0 under the existing +repository notice. No license change or exclusive rights are asserted. + +This preview uses NumPy, Matplotlib and tmm as separately installed +dependencies. Preserve their respective notices when redistributing them. +The preview wheel does not bundle those dependencies. + +Method reference: +S. J. Byrnes, "Multilayer optical calculations", arXiv:1603.02720. +https://doi.org/10.48550/arXiv.1603.02720 + +Material index formulas are inherited illustrative placeholders with no +measurement provenance. This software is not a validated OLED dipole +emission/outcoupling solver or a certified device-design tool. diff --git a/README.md b/README.md index 8689486..3f01f16 100644 --- a/README.md +++ b/README.md @@ -2,13 +2,129 @@ Python scripts and Mathematica notebook revisions exploring multilayer optics, polarization, wavelength dependence, and OLED outcoupling concepts. +## Installable plane-wave preview (0.1.0a1) + +**Release preparation only; not approved for paid publication.** The packaged +program is a headless plane-wave transmittance demonstration, not an OLED +dipole-emission, outcoupling-efficiency, chromaticity, or device-design solver. +The inherited organic and silver index formulas are illustrative placeholders, +not measured material data. Historical scripts and notebooks are not included +in the wheel and are not certified by these changes. + +Python 3.11-3.14 is the declared compatibility range; CI checks the matrix below, +and only completed test runs establish support. From the repository, on Windows: + +```powershell +python -m venv .venv +.\.venv\Scripts\python.exe -m pip install ".[dev]" +.\.venv\Scripts\python.exe -m pytest -q +.\.venv\Scripts\mongelli-plane-wave.exe --wavelength-count 31 --angle-count 19 --output transmission.png +``` + +The command creates: + +- `transmission.png`: a headless wavelength/angle contour plot. +- `transmission.csv`: columns `wavelength_nm`, `angle_deg`, + `unpolarized_transmittance`; wavelength-major order, one row per sample. +- `transmission.json`: software version, coordinates, layer thicknesses, + polarization normalization, and model limitations. + +Lengths use nm, angles use degrees, and transmittance is dimensionless. +The stack is air / illustrative organic / 50 nm illustrative silver / air, +with semi-infinite air on both sides. The default grid is 400-700 nm and +0-89 degrees; unpolarized transmittance is `(T_s + T_p) / 2`, not their sum. +Use `--thickness-nm` to change the organic thickness. Counts must be at least +two. Output must end in `.png` and its parent directory must already exist. +Existing output files are protected unless `--overwrite` is explicitly passed. +Errors exit nonzero; a write failure may leave partial outputs. + +For a built wheel, install with `python -m pip install ` and use +`mongelli-plane-wave` from the same environment. An internet connection is +needed to resolve dependencies unless wheels were separately supplied. This +is **not** a standalone Windows installer, and does not bundle Mathematica, +NumPy, Matplotlib, or tmm. + +### Validation and distribution + +```powershell +.\.venv\Scripts\python.exe -m build +.\.venv\Scripts\python.exe -m twine check dist\* +``` + +Tests cover analytic Fresnel interfaces, Airy thin-film transmittance, a separate +normal-incidence characteristic-matrix calculation for the absorbing stack, +lossless power balance, invalid inputs, and CLI export. These checks validate +bounded plane-wave behavior; they do not validate emitter physics or measured +material properties. See `tests/test_tmm_oled9.py` for cases and tolerances. + +The GitHub workflow runs numerical tests, builds distributions, smoke-tests an +installed wheel outside the checkout, and uploads preview artifacts. Its +runtime-dependency audit emits a CycloneDX inventory and fails on known +vulnerabilities (or unavailable advisory data); it is not a full security review. +The inventory also includes environment tooling such as pip. Its +optional manual draft-release action accepts only a tag matching the package +version and creates a **draft prerelease**; it never publishes a paid listing. +No automatic tag push or public release is performed. + +### GitHub + Gumroad commercial approval gates + +The five enrollment packets are supporting checklists, not enrollment approvals. +Use GitHub for technical documentation, demonstrations and versioned artifacts; +Gumroad is the proposed separate checkout and buyer-delivery channel. + +Before charging customers: + +1. Decide whether the product is this educational plane-wave tool or a future + validated OLED solver. Do not use OLED efficiency claims for this preview. +2. Confirm code/contributor and material-data provenance and preserve the + existing Apache-2.0 code and CC BY 4.0 documentation terms. A purchase cannot + remove existing open-source rights; exclusive licensing is not implied. +3. Complete independent scientific review for the intended claims, dependency + vulnerability scanning, license review, supported-platform testing, release + checksums, a dependency SBOM, and signed-tag provenance. +4. Approve price, customer license wording, support contact and response policy, + refund terms, updates policy, privacy notice, and truthful screenshots. +5. Complete seller identity, tax and payout verification directly in Gumroad. + Never store tax identifiers, identity documents, or banking data here. +6. Test checkout and buyer delivery using the platform's supported test process. + Upload paid deliverables to authorized buyer delivery, not a public GitHub + release if they are intended to be access-restricted. + +No price, exclusive license, seller approval, support SLA, paid listing, or +production readiness is asserted by this repository. Microsoft/AWS/GCP billing +integrations are intentionally outside this first downloadable-preview scope. + +### Proposed storefront copy (draft, not a live offer) + +**Title:** Mongelli Plane-Wave Optics Preview + +**Audience:** Researchers, educators and scientific programmers exploring +planar thin-film transmission with Python. + +**Description:** Generate wavelength/angle transmission maps for an illustrative +air/organic/silver/air stack. Export a PNG plot, numerical CSV and model metadata. +Includes source code and a Python command-line interface. Requires Python and +separately installed dependencies; no standalone Windows executable is supplied. + +**Limitations:** Alpha research/educational preview. No measured material +database, emitting dipoles, calibrated OLED efficiency, GUI, or device-performance +guarantee. The independent checks cover specified plane-wave cases only. + +**Rights:** Existing Apache-2.0 code and CC BY 4.0 documentation terms remain. +Any eventual charge must accurately describe the added value (for example, +packaging or a separately agreed support service), not imply exclusive rights. + +**Not yet approved:** Price, support commitments, refund/updates terms, customer +license presentation, final screenshots, seller verification and publication. + ## Functions and engineering applications | Source / interface | Observed function | Inferred application | |---|---|---| | `python_tmm_oled_Green'sfcnApproach.py`: `OLEDOutcouplingTMM` | Layer definitions, aluminum Drude model, interface/propagation matrices, angle/wavelength calculations | Exploring material and thickness sensitivity in OLED stacks [1] | | `oled_outcoupling_simulation.py`: `transfer_matrix_method(...)` | Normal-incidence characteristic-matrix calculation | Teaching and checking thin-film interference concepts [1] | -| `tmm_OLED8.py`, `tmm_OLED9.py` | Wavelength/angle/thickness sweeps and image output | Exploratory design maps, after correcting implementation issues | +| `tmm_OLED9.py` | Installable plane-wave transmittance sweep, PNG/CSV/JSON export | Illustrative optics only; not OLED dipole outcoupling | +| `tmm_OLED8.py` | Historical wavelength/angle/thickness sweep | Unvalidated prototype, excluded from the package | | `TMMandGreensFunctions*.nb`, `mathematica_tmm_oled.txt` | Symbolic/numerical multilayer and Green-function formulations | Comparing analytical formulations with numerical solvers [1,2] | These are research prototypes. The class filename does not establish a complete dyadic Green-function emission calculation: its efficiency method explicitly uses a simplified transmission-amplitude approach. diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..3794ce0 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,39 @@ +[build-system] +requires = ["setuptools>=77"] +build-backend = "setuptools.build_meta" + +[project] +name = "mongelli-plane-wave-preview" +version = "0.1.0a1" +description = "Illustrative multilayer plane-wave optics preview, not an OLED emission solver" +readme = "README.md" +requires-python = ">=3.11,<3.15" +license = "Apache-2.0" +license-files = ["LICENSE", "LICENSE-CC-BY-4.0.txt", "NOTICE"] +authors = [{name = "Guy Francis Mongelli"}] +dependencies = [ + "numpy>=1.24,<3", + "matplotlib>=3.7,<4", + "tmm>=0.2.0,<0.3", +] +classifiers = [ + "Development Status :: 3 - Alpha", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Physics", +] + +[project.optional-dependencies] +dev = ["pytest>=8,<10", "build>=1.2,<2", "twine>=6,<7", "pip-audit>=2.9,<3"] + +[project.scripts] +mongelli-plane-wave = "tmm_OLED9:main" + +[project.urls] +Repository = "https://github.com/gmongell/Mathematica_python_TMM_OLEDPixel" + +[tool.setuptools] +py-modules = ["tmm_OLED9"] + +[tool.pytest.ini_options] +testpaths = ["tests"] +pythonpath = ["."] diff --git a/tests/test_tmm_oled9.py b/tests/test_tmm_oled9.py new file mode 100644 index 0000000..a0770dc --- /dev/null +++ b/tests/test_tmm_oled9.py @@ -0,0 +1,213 @@ +import csv +import json +import os +from pathlib import Path +import subprocess +import sys + +import numpy as np +import pytest + +import tmm_OLED9 as model + + +@pytest.mark.parametrize("polarization", ["s", "p"]) +@pytest.mark.parametrize("angle", [0.0, 30.0, 65.0, 89.0]) +def test_zero_thickness_is_air(polarization, angle): + result = model.power_coefficients(550.0, angle, 0.0, polarization, silver_nm=0.0) + assert result["R"] == pytest.approx(0.0, abs=1e-12) + assert result["T"] == pytest.approx(1.0, abs=1e-12) + assert result["A"] == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("polarization", ["s", "p"]) +@pytest.mark.parametrize("angle", [0.0, 30.0, 65.0]) +def test_interface_matches_fresnel(monkeypatch, polarization, angle): + monkeypatch.setattr(model, "stack_at", lambda *args: ([1.0, 1.5], [np.inf, np.inf])) + theta = np.deg2rad(angle) + incident_cos = np.cos(theta) + transmitted_cos = np.sqrt(1 - (np.sin(theta) / 1.5) ** 2) + if polarization == "s": + r = (incident_cos - 1.5 * transmitted_cos) / (incident_cos + 1.5 * transmitted_cos) + else: + r = (1.5 * incident_cos - transmitted_cos) / (1.5 * incident_cos + transmitted_cos) + result = model.power_coefficients(550.0, angle, 0.0, polarization) + assert result["R"] == pytest.approx(r**2, abs=1e-12) + assert result["T"] == pytest.approx(1 - r**2, abs=1e-12) + + +@pytest.mark.parametrize("polarization", ["s", "p"]) +def test_lossless_film_matches_airy_formula(polarization): + wavelength, thickness = 550.0, 75.0 + index = 1.65 + 0.5 / 550.0**2 + interface_reflectance = ((index - 1) / (index + 1)) ** 2 + phase = 2 * np.pi * index * thickness / wavelength + expected = 1 / ( + 1 + 4 * interface_reflectance * np.sin(phase) ** 2 / (1 - interface_reflectance) ** 2 + ) + result = model.power_coefficients(wavelength, 0.0, thickness, polarization, silver_nm=0) + assert result["T"] == pytest.approx(expected, abs=1e-12) + assert result["R"] == pytest.approx(1 - expected, abs=1e-12) + assert result["A"] == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("polarization", ["s", "p"]) +def test_absorbing_stack_matches_characteristic_matrix(polarization): + # Independent normal-incidence field-matrix formulation, not another coh_tmm call. + matrix = np.eye(2, dtype=complex) + for index, thickness in [(1.65 + 0.5 / 550**2, 75), (0.16 - 1.5 / 550**2 + 3.16j, 50)]: + phase = 2 * np.pi * index * thickness / 550 + matrix = matrix @ np.array([ + [np.cos(phase), -1j * np.sin(phase) / index], + [-1j * index * np.sin(phase), np.cos(phase)], + ]) + b = matrix[0, 0] + matrix[0, 1] + c = matrix[1, 0] + matrix[1, 1] + expected_r = abs((b - c) / (b + c)) ** 2 + expected_t = abs(2 / (b + c)) ** 2 + result = model.power_coefficients(550.0, 0.0, 75.0, polarization) + assert result["R"] == pytest.approx(expected_r, abs=1e-12) + assert result["T"] == pytest.approx(expected_t, abs=1e-12) + assert result["A"] == pytest.approx(1 - expected_r - expected_t, abs=1e-12) + assert 0 < result["A"] < 1 + + +@pytest.mark.parametrize("polarization", ["s", "p"]) +def test_lossless_oblique_power_balance(polarization): + for wavelength in [400.0, 550.0, 700.0]: + for angle in [0.0, 45.0, 89.0]: + result = model.power_coefficients(wavelength, angle, 100.0, polarization, silver_nm=0) + assert result["R"] + result["T"] == pytest.approx(1.0, abs=1e-12) + assert result["A"] == pytest.approx(0.0, abs=1e-12) + + +def test_sweep_shape_weighting_and_high_angle_transmission(): + values = model.sweep([450.0, 550.0], [0.0, 65.0], 75.0) + assert values.shape == (2, 2) + for i, wavelength in enumerate([450.0, 550.0]): + for j, angle in enumerate([0.0, 65.0]): + expected = sum(model.transmission(wavelength, angle, 75.0, p) for p in ("s", "p")) / 2 + assert values[i, j] == pytest.approx(expected, abs=1e-12) + assert np.all(values > 0) + assert np.all(values <= 1) + + +@pytest.mark.parametrize("wavelength,angle,organic,silver,polarization", [ + (0, 0, 75, 50, "s"), (-1, 0, 75, 50, "s"), (np.nan, 0, 75, 50, "s"), + (np.inf, 0, 75, 50, "s"), (550, -1, 75, 50, "s"), (550, 90, 75, 50, "s"), + (550, np.nan, 75, 50, "s"), (550, np.inf, 75, 50, "s"), + (550, 0, -1, 50, "s"), (550, 0, np.nan, 50, "s"), + (550, 0, np.inf, 50, "s"), (550, 0, 75, -1, "s"), + (550, 0, 75, np.nan, "s"), (550, 0, 75, np.inf, "s"), (550, 0, 75, 50, "x"), +]) +def test_invalid_scalar_inputs(wavelength, angle, organic, silver, polarization): + with pytest.raises(ValueError): + model.power_coefficients(wavelength, angle, organic, polarization, silver) + + +@pytest.mark.parametrize("wavelengths,angles,thickness", [ + ([], [0], 75), ([550], [], 75), ([[550]], [0], 75), ([550], [[0]], 75), + ([0], [0], 75), ([np.nan], [0], 75), ([550], [90], 75), + ([550], [np.inf], 75), ([550], [0], -1), +]) +def test_invalid_sweep_inputs(wavelengths, angles, thickness): + with pytest.raises(ValueError): + model.sweep(wavelengths, angles, thickness) + + +def test_cli_exports_plot_csv_and_metadata(tmp_path, capsys): + output = tmp_path / "preview.png" + model.main(["--wavelength-count", "2", "--angle-count", "3", "--output", str(output)]) + assert output.read_bytes().startswith(b"\x89PNG\r\n\x1a\n") + with output.with_suffix(".csv").open(newline="", encoding="utf-8") as stream: + rows = list(csv.DictReader(stream)) + assert len(rows) == 6 + assert [(float(row["wavelength_nm"]), float(row["angle_deg"])) for row in rows] == [ + (400, 0), (400, 44.5), (400, 89), (700, 0), (700, 44.5), (700, 89), + ] + expected = model.sweep([400, 700], [0, 44.5, 89], 75) + np.testing.assert_allclose( + [float(row["unpolarized_transmittance"]) for row in rows], expected.ravel(), atol=1e-12 + ) + metadata = json.loads(output.with_suffix(".json").read_text(encoding="utf-8")) + assert metadata["quantity"] == "unpolarized_plane_wave_transmittance" + assert metadata["organic_thickness_nm"] == 75 + assert metadata["silver_thickness_nm"] == 50 + assert metadata["wavelengths_nm"] == [400, 700] + assert metadata["angles_deg"] == [0, 44.5, 89] + assert metadata["software_version"] == model.__version__ + assert "NOT a validated OLED solver" in capsys.readouterr().err + + +@pytest.mark.parametrize("arguments,message", [ + (["--wavelength-count", "1"], "sample counts"), + (["--angle-count", "1"], "sample counts"), + (["--thickness-nm", "nan"], "finite and nonnegative"), + (["--thickness-nm", "-1"], "finite and nonnegative"), +]) +def test_cli_rejects_invalid_inputs(tmp_path, capsys, arguments, message): + with pytest.raises(SystemExit) as error: + model.main(arguments + ["--output", str(tmp_path / "preview.png")]) + assert error.value.code != 0 + assert message in capsys.readouterr().err + assert list(tmp_path.iterdir()) == [] + + +def test_cli_does_not_overwrite_existing_sidecar(tmp_path, capsys): + output = tmp_path / "preview.png" + output.with_suffix(".csv").write_text("keep me", encoding="utf-8") + with pytest.raises(SystemExit) as error: + model.main(["--output", str(output)]) + assert error.value.code == 2 + assert "output already exists" in capsys.readouterr().err + assert output.with_suffix(".csv").read_text(encoding="utf-8") == "keep me" + assert not output.exists() + + +def test_cli_explicit_overwrite_replaces_outputs(tmp_path): + output = tmp_path / "preview.png" + for suffix in [".png", ".csv", ".json"]: + output.with_suffix(suffix).write_text("old content", encoding="utf-8") + model.main([ + "--wavelength-count", "2", "--angle-count", "2", + "--output", str(output), "--overwrite", + ]) + assert output.read_bytes().startswith(b"\x89PNG\r\n\x1a\n") + assert output.with_suffix(".csv").read_text(encoding="utf-8").startswith("wavelength_nm,") + assert json.loads(output.with_suffix(".json").read_text(encoding="utf-8"))["schema_version"] == 1 + + +@pytest.mark.parametrize("relative_path,message", [ + ("bad.svg", ".png extension"), + ("missing/preview.png", "parent directory"), +]) +def test_cli_rejects_invalid_output_paths(tmp_path, capsys, relative_path, message): + with pytest.raises(SystemExit) as error: + model.main(["--output", str(tmp_path / relative_path)]) + assert error.value.code != 0 + assert message in capsys.readouterr().err + assert list(tmp_path.iterdir()) == [] + + +def test_cli_reports_write_failure(tmp_path, monkeypatch, capsys): + def fail_open(*args, **kwargs): + raise PermissionError("test write denied") + + monkeypatch.setattr(Path, "open", fail_open) + with pytest.raises(SystemExit) as error: + model.main([ + "--wavelength-count", "2", "--angle-count", "2", + "--output", str(tmp_path / "preview.png"), + ]) + assert error.value.code == 1 + assert "test write denied" in capsys.readouterr().err + + +def test_import_has_no_plot_or_file_side_effects(tmp_path): + env = dict(os.environ, PYTHONPATH=str(Path(__file__).resolve().parents[1])) + result = subprocess.run( + [sys.executable, "-c", "import sys; import tmm_OLED9; assert 'matplotlib.pyplot' not in sys.modules"], + cwd=tmp_path, env=env, capture_output=True, text=True, check=False, + ) + assert result.returncode == 0, result.stderr + assert list(tmp_path.iterdir()) == [] diff --git a/tmm_OLED9.py b/tmm_OLED9.py index 3208297..59e7364 100644 --- a/tmm_OLED9.py +++ b/tmm_OLED9.py @@ -1,80 +1,174 @@ -import numpy as np -import matplotlib -matplotlib.use('Agg') # Use Agg backend for headless environments -import matplotlib.pyplot as plt - -# Dispersion relationships for materials -def n_alq3(wavelength): - A = 1.65 - B = 0.5 # Adjust as necessary - return A + B / (wavelength ** 2) - -def n_silver(wavelength): - n_0 = 0.16 # Real part - kappa = 3.16 # Imaginary part - C = 1.5 # Coefficient for the dispersion - return n_0 - (C / (wavelength ** 2)) + 1j * kappa - -# Wavelength and angle ranges -wavelengths = np.linspace(400, 700, 1000) # Wavelength range (in nm) -theta_list_deg = np.linspace(0, 90, 180) # Angles from 0 to 90 degrees -theta_list_rad = np.radians(theta_list_deg) # Convert to radians - -# Set thickness range for the organic layer (e.g., Alq3) in nanometers -thicknesses = np.linspace(25, 150, 10) # Example thickness values, from 50nm to 200nm - -# Loop over different thickness values of the organic layer -for thickness in thicknesses: - # Initialize the outcoupling efficiency array for this thickness - outcoupling_efficiency_2D = np.zeros((len(wavelengths), len(theta_list_deg))) - - # Create the n_list for the current wavelength range and thickness - n_air = 1.0 # Refractive index of air - n_list = [n_air] # Start with air - for wavelength in wavelengths: - n_list.append(n_alq3(wavelength)) # Alq3 layer refractive index - n_list.append(n_silver(wavelengths[-1])) # Silver refractive index at last wavelength - - # Set the layer thickness values - d_alq3 = thickness * 1e-9 # Convert nm to meters - d_silver = 50e-9 # 50 nm for Silver - d_list = [np.inf, d_alq3, d_silver, np.inf] # Air, Alq3, Silver, and back air - - # Loop over each wavelength and angle to compute the outcoupling efficiency - for i, wavelength in enumerate(wavelengths): - R_TE, T_TE, R_TM, T_TM = [], [], [], [] - - for theta in theta_list_rad: - # Ensure angle is valid for first layer (Alq3) - if theta < np.arcsin(n_air / np.abs(n_list[1])): - # TE polarization (s-polarized light) - coh_TE = coh_tmm('s', n_list, d_list, wavelength * 1e-9, theta) - R_TE.append(coh_TE['R']) - T_TE.append(coh_TE['T']) - - # TM polarization (p-polarized light) - coh_TM = coh_tmm('p', n_list, d_list, wavelength * 1e-9, theta) - R_TM.append(coh_TM['R']) - T_TM.append(coh_TM['T']) - else: - R_TE.append(1.0) # Total reflection - T_TE.append(0.0) # No transmission - R_TM.append(1.0) # Total reflection - T_TM.append(0.0) # No transmission - - # Combine TE and TM transmission efficiencies to compute the outcoupling efficiency - outcoupling_efficiency = np.array(T_TE) + np.array(T_TM) - outcoupling_efficiency_2D[i, :] = outcoupling_efficiency - - # Plotting the outcoupling efficiency for this thickness as a heatmap - plt.figure(figsize=(10, 6)) - plt.contourf(theta_list_deg, wavelengths, outcoupling_efficiency_2D, levels=50, cmap='viridis') - plt.colorbar(label='Outcoupling Efficiency') - plt.xlabel('Angle of Incidence (degrees)') - plt.ylabel('Wavelength (nm)') - plt.title(f'Outcoupling Efficiency (Alq3 Thickness = {thickness:.1f} nm)') - - # Save the plot as an image file - plt.savefig(f'outcoupling_efficiency_thickness_{thickness:.1f}nm.png') - plt.close() # Close the figure to free up memory - +"""Illustrative plane-wave multilayer optics; NOT OLED dipole outcoupling.""" + +import argparse +import csv +import json +import sys +from pathlib import Path +from typing import Sequence + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from tmm import coh_tmm + +__version__ = "0.1.0a1" + + +def _positive_wavelength(wavelength_nm: float) -> None: + if not np.isfinite(wavelength_nm) or wavelength_nm <= 0: + raise ValueError("Wavelength must be finite and positive (nm)") + + +def _nonnegative_thickness(thickness_nm: float) -> None: + if not np.isfinite(thickness_nm) or thickness_nm < 0: + raise ValueError("Thickness must be finite and nonnegative (nm)") + + +def n_alq3(wavelength_nm: float) -> float: + """Illustrative dispersion only, not measured or calibrated Alq3 data.""" + _positive_wavelength(wavelength_nm) + return 1.65 + 0.5 / wavelength_nm**2 + + +def n_silver(wavelength_nm: float) -> complex: + """Illustrative passive index, not measured or calibrated silver data.""" + _positive_wavelength(wavelength_nm) + return 0.16 - 1.5 / wavelength_nm**2 + 3.16j + + +def stack_at( + wavelength_nm: float, organic_nm: float, silver_nm: float = 50.0 +) -> tuple[list[complex], list[float]]: + _positive_wavelength(wavelength_nm) + _nonnegative_thickness(organic_nm) + _nonnegative_thickness(silver_nm) + return ( + [1.0, n_alq3(wavelength_nm), n_silver(wavelength_nm), 1.0], + [np.inf, organic_nm, silver_nm, np.inf], + ) + + +def power_coefficients( + wavelength_nm: float, + angle_deg: float, + organic_nm: float, + polarization: str, + silver_nm: float = 50.0, +) -> dict[str, float]: + """Return R, T and total A for air/organic/silver/air plane-wave incidence. + + A = 1 - R - T applies here because both external media are lossless air. + All lengths are nm; angle is in degrees. No emitting dipoles are modeled. + """ + if not np.isfinite(angle_deg) or not 0 <= angle_deg < 90: + raise ValueError("Angle must be finite and in [0, 90) degrees") + if polarization not in ("s", "p"): + raise ValueError("Polarization must be s or p") + indices, thicknesses = stack_at(wavelength_nm, organic_nm, silver_nm) + result = coh_tmm( + polarization, indices, thicknesses, np.deg2rad(angle_deg), wavelength_nm + ) + reflection, transmittance = float(result["R"]), float(result["T"]) + absorption = 1.0 - reflection - transmittance + values = np.array([reflection, transmittance, absorption]) + if not np.all(np.isfinite(values)) or np.any(values < -1e-9) or np.any(values > 1 + 1e-9): + raise ArithmeticError("Solver returned nonphysical passive-stack power coefficients") + return {"R": reflection, "T": transmittance, "A": absorption} + + +def transmission( + wavelength_nm: float, angle_deg: float, organic_nm: float, polarization: str +) -> float: + return power_coefficients(wavelength_nm, angle_deg, organic_nm, polarization)["T"] + + +def sweep( + wavelengths_nm: ArrayLike, angles_deg: ArrayLike, organic_nm: float +) -> NDArray[np.float64]: + """Return equally weighted s/p transmittance, indexed [wavelength, angle].""" + wavelengths = np.asarray(wavelengths_nm, dtype=float) + angles = np.asarray(angles_deg, dtype=float) + _nonnegative_thickness(organic_nm) + if wavelengths.ndim != 1 or angles.ndim != 1 or not wavelengths.size or not angles.size: + raise ValueError("Sweep coordinates must be nonempty 1-D arrays") + if not np.all(np.isfinite(wavelengths)) or np.any(wavelengths <= 0): + raise ValueError("Wavelengths must be finite and positive (nm)") + if not np.all(np.isfinite(angles)) or np.any(angles < 0) or np.any(angles >= 90): + raise ValueError("Angles must be finite and in [0, 90) degrees") + values = np.empty((len(wavelengths), len(angles))) + for i, wavelength in enumerate(wavelengths): + for j, angle in enumerate(angles): + values[i, j] = 0.5 * ( + transmission(wavelength, angle, organic_nm, "s") + + transmission(wavelength, angle, organic_nm, "p") + ) + return values + + +def main(argv: Sequence[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--version", action="version", version=__version__) + parser.add_argument("--wavelength-count", type=int, default=31) + parser.add_argument("--angle-count", type=int, default=19) + parser.add_argument("--thickness-nm", type=float, default=75.0) + parser.add_argument("--output", type=Path, default=Path("transmission.png")) + parser.add_argument("--overwrite", action="store_true", help="Replace existing output files") + args = parser.parse_args(argv) + if args.wavelength_count < 2 or args.angle_count < 2: + parser.error("sample counts must be >= 2") + if args.output.suffix.lower() != ".png": + parser.error("--output must have a .png extension") + paths = [args.output, args.output.with_suffix(".csv"), args.output.with_suffix(".json")] + if not args.overwrite and any(path.exists() for path in paths): + parser.error("output already exists; choose another path or use --overwrite") + try: + _nonnegative_thickness(args.thickness_nm) + if not args.output.parent.is_dir(): + raise ValueError("Output parent directory must already exist") + wavelengths = np.linspace(400.0, 700.0, args.wavelength_count) + angles = np.linspace(0.0, 89.0, args.angle_count) + values = sweep(wavelengths, angles, args.thickness_nm) + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, ax = plt.subplots() + try: + plot = ax.contourf(angles, wavelengths, values, levels=30) + fig.colorbar(plot, ax=ax, label="Unpolarized plane-wave transmittance") + ax.set( + xlabel="Incidence angle (degrees)", + ylabel="Vacuum wavelength (nm)", + title=f"Illustrative stack: organic {args.thickness_nm:g} nm", + ) + fig.savefig(paths[0], dpi=150) + finally: + plt.close(fig) + with paths[1].open("w", newline="", encoding="utf-8") as output: + writer = csv.writer(output) + writer.writerow(["wavelength_nm", "angle_deg", "unpolarized_transmittance"]) + for i, wavelength in enumerate(wavelengths): + for j, angle in enumerate(angles): + writer.writerow([wavelength, angle, values[i, j]]) + metadata = { + "schema_version": 1, + "software_version": __version__, + "quantity": "unpolarized_plane_wave_transmittance", + "polarization_weighting": "0.5 * (T_s + T_p)", + "stack": ["air", "illustrative organic", "illustrative silver", "air"], + "organic_thickness_nm": args.thickness_nm, + "silver_thickness_nm": 50.0, + "wavelengths_nm": wavelengths.tolist(), + "angles_deg": angles.tolist(), + "limitations": "Illustrative indices; not calibrated materials or OLED dipole outcoupling.", + } + paths[2].write_text(json.dumps(metadata, indent=2) + "\n", encoding="utf-8") + except (ValueError, ArithmeticError, OSError) as error: + parser.exit(1, f"error: {error}; outputs may be incomplete\n") + print("WARNING: illustrative plane-wave preview; NOT a validated OLED solver.", file=sys.stderr) + print("Created: " + ", ".join(str(path) for path in paths)) + + +if __name__ == "__main__": + main() From 58c8ab3dd7aff213e34886d018d828a6dfe0dce6 Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 17:06:24 -0400 Subject: [PATCH 2/7] ci: refresh buyer installer and retain dependency audit evidence Upgrade the venv bootstrap pip before clean installation and retain the CycloneDX report even when auditing blocks distribution. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .github/workflows/plane-wave-preview.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/plane-wave-preview.yml b/.github/workflows/plane-wave-preview.yml index 53d92cd..bd75a2a 100644 --- a/.github/workflows/plane-wave-preview.yml +++ b/.github/workflows/plane-wave-preview.yml @@ -55,6 +55,7 @@ jobs: shell: bash run: | python -m venv "$RUNNER_TEMP/buyer" + "$RUNNER_TEMP/buyer/bin/python" -m pip install --upgrade pip "$RUNNER_TEMP/buyer/bin/python" -m pip install dist/*.whl cd "$RUNNER_TEMP" "$RUNNER_TEMP/buyer/bin/mongelli-plane-wave" --wavelength-count 3 --angle-count 3 --output preview.png @@ -63,6 +64,12 @@ jobs: run: | SITE=$("$RUNNER_TEMP/buyer/bin/python" -c "import sysconfig; print(sysconfig.get_paths()['purelib'])") python -m pip_audit --path "$SITE" --format cyclonedx-json --output dist/dependency-audit.cdx.json + - uses: actions/upload-artifact@v4 + if: always() + with: + name: runtime-dependency-audit + path: dist/dependency-audit.cdx.json + if-no-files-found: warn - name: Record artifact checksums run: sha256sum dist/*.whl dist/*.tar.gz dist/*.json > dist/SHA256SUMS.txt - uses: actions/upload-artifact@v4 From 0e5bfd7c5a94fe51ef8b8aa68f525b1f4ba6f471 Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 17:24:10 -0400 Subject: [PATCH 3/7] docs: draft seller policies and OLED solver acceptance scope Prepare owner-review pricing, support/refund/privacy copy, provenance gates and verified Gumroad onboarding guidance. Record Edge login handoff and specify dipole solver physics and acceptance tests without claiming implementation or seller approval. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- README.md | 6 + docs/COMMERCIAL_DRAFTS.md | 212 +++++++++++++++++++++++++++++++ docs/GUMROAD_SELLER_CHECKLIST.md | 159 +++++++++++++++++++++++ docs/OLED_SOLVER_SPEC.md | 151 ++++++++++++++++++++++ 4 files changed, 528 insertions(+) create mode 100644 docs/COMMERCIAL_DRAFTS.md create mode 100644 docs/GUMROAD_SELLER_CHECKLIST.md create mode 100644 docs/OLED_SOLVER_SPEC.md diff --git a/README.md b/README.md index 3f01f16..6858c37 100644 --- a/README.md +++ b/README.md @@ -68,6 +68,12 @@ No automatic tag push or public release is performed. ### GitHub + Gumroad commercial approval gates +Owner-review drafts are available in [commercial policies and storefront copy](docs/COMMERCIAL_DRAFTS.md), +[Gumroad seller verification guidance](docs/GUMROAD_SELLER_CHECKLIST.md), and +[the proposed OLED dipole-solver specification](docs/OLED_SOLVER_SPEC.md). +These documents are not active policies, completed verification, or implemented +solver functionality. + The five enrollment packets are supporting checklists, not enrollment approvals. Use GitHub for technical documentation, demonstrations and versioned artifacts; Gumroad is the proposed separate checkout and buyer-delivery channel. diff --git a/docs/COMMERCIAL_DRAFTS.md b/docs/COMMERCIAL_DRAFTS.md new file mode 100644 index 0000000..7c256a1 --- /dev/null +++ b/docs/COMMERCIAL_DRAFTS.md @@ -0,0 +1,212 @@ +# Commercial launch drafts + +Prepared 2026-10-16. **DRAFT: not approved, published, or binding.** +These proposed terms need owner review and appropriate legal/tax advice. +Do not copy unresolved fields into a live listing. No seller registration, +professional qualification, scientific certification, or sales is asserted. + +## 1. Product and pricing proposal + +Proposed public brand: Mongelli Labs, subject to trade-name verification. +The contracting seller's verified legal identity must be supplied through +the appropriate platform fields before publication. + +| Item | Proposal | Release condition | +|---|---|---| +| Current alpha | Free technical evaluation; no paid launch now | Owner approves any public release; prominently retain limitations | +| Later educational bundle | USD 19 one-time starting-price hypothesis | Scientific scope, provenance, installation and commercial gates pass | +| Future OLED solver | Not priced or offered yet | Separate implementation, validation, scope and support review | +| Consulting/custom development | Separate written quote | No engineering certification, outcome or delivery date implied | + +USD 19 is a positioning experiment, not market-researched willingness to pay. +Charge for a documented, convenient bundle or expressly agreed services, +not exclusive access to existing Apache-2.0 source. Do not imply a purchase +is required to exercise rights already granted by the open-source licenses. +Do not promise the future OLED solver as part of this purchase. + +### Fee illustration, not net-income advice + +Official Gumroad fees checked on 2026-10-16: + +- Direct website sale: 10% + USD 0.50, with credit-card processing separately + stated as 2.9% + USD 0.30. PayPal fees vary. +- Discover marketplace sale: flat 30%, including processing. +- Different account-specific rates, applicable taxes, currency/payout charges, + refunds, and future fee changes can alter proceeds. + +For a USD 19 direct card sale, the illustrative remainder is +`19 - (19 * 0.10 + 0.50) - (19 * 0.029 + 0.30) = USD 15.749`, +approximately USD 15.75 before other applicable adjustments and business taxes. +At 30% Discover fees, the comparable amount is USD 13.30. +Neither figure is a guaranteed payout or profit. Recheck live fees at launch. + +Gumroad's merchant-of-record sales-tax handling does not eliminate the seller's +income/self-employment taxes, recordkeeping or other business obligations. + +## 2. Storefront copy + +**Title:** Mongelli Plane-Wave Optics Educational Preview + +**Summary:** Explore thin-film transmission with a Python command-line tool. +Export numerical results, a contour plot, and reproducibility metadata. + +**For:** Researchers, educators, and scientific programmers learning planar +multilayer optics. This is not a design-certified device simulator. + +**Included in the proposed bundle:** + +- Versioned Python wheel and source archive with license notices. +- Installation guide and command-line examples. +- Example PNG, CSV and JSON outputs. +- Validation notes covering specified plane-wave reference cases. +- Release checksums and runtime-dependency CycloneDX inventory. + +**Requirements:** Python 3.11-3.14, ability to install Python packages, and +internet access to resolve dependencies unless separately supplied. Tests have +passed on GitHub-hosted Windows and Linux across that Python range. +macOS has not been tested. No standalone Windows installer or GUI is supplied. + +**Model:** Air / illustrative organic film / 50 nm illustrative silver / air. +Default wavelengths: 400-700 nm; incidence angles: 0-89 degrees. +Transmittance is the equal-weight average of s and p polarizations. + +**Limitations:** Alpha research/educational software. The index formulas are +illustrative and have no verified measurement provenance. No emitting dipoles, +calibrated OLED outcoupling, external quantum efficiency, device lifetime, +color accuracy, or fabrication-performance guarantee. Historical notebooks +are not part of the installable wheel. Do not describe the preview as the +future OLED solver. + +**Price/support/refunds:** Insert only the owner-approved price and policies. +Keep this product unpublished until every required gate is complete. + +## 3. Proposed customer terms + +### Software rights and intended use + +Project code remains Apache-2.0; project documentation and figures remain +CC BY 4.0 under the existing repository notices. Third-party components retain +their respective licenses. Preserve applicable attribution and notices. +The purchase, if offered, is nonexclusive and does not revoke or replace those +licenses. The educational description states the intended scope; it does not +add a field-of-use restriction to Apache-2.0. + +No independently certified engineering result or device-performance warranty +is offered. Existing license disclaimers apply to the extent permitted by law; +mandatory consumer rights are not excluded. No additional bespoke liability +cap, arbitration clause, or governing-law clause is proposed here. + +### Proposed support and updates policy + +- Private support: use the seller contact shown on the Gumroad receipt. + Confirm the mailbox and its monitoring before publishing. +- Public, non-sensitive reproducible bugs: use the repository issue tracker. +- Proposed response target: within three business days, Monday-Friday, + excluding seller-observed US holidays. This is a target, not a guaranteed SLA. +- Proposed installation assistance window: 30 calendar days from purchase, + limited to documented commands on supported platforms. +- No remote machine administration, confidential-data analysis, custom modeling, + or scientific certification is included. +- Purchase delivers the stated version. Future fixes may be released, but + perpetual support, all future versions, or the OLED solver are not promised. +- Announce critical known defects through available release/advisory channels; + do not promise uninterrupted service or defect-free results. + +These are proposed commitments only. The owner must confirm support capacity +before activation, including arrangements for absences or discontinuation. + +### Proposed refund policy + +Suggested customer-facing text, effective only after owner approval: + +> Contact the seller through your Gumroad receipt within 30 calendar days of +> purchase to request a refund if the download cannot be used on a documented +> supported configuration or the product materially differs from its listing. +> Please provide your receipt reference, operating system, Python version and +> a non-sensitive error description. Troubleshooting is optional; we will not +> require confidential files. Refund requests are handled through Gumroad. +> Mandatory statutory rights and Gumroad's applicable policies remain in effect, +> including any legally required withdrawal rights. + +Do not require deletion of Apache-licensed code or surrender of existing +license rights as a refund condition. Refund timing and payment-processing +charges are governed by the platform and payment rail, not guaranteed here. +Confirm the exact policy in the product's refund settings before launch. + +### Proposed privacy notice + +Draft operational commitments, not a statement that a live service exists: + +> This product runs locally. The packaged preview has no application-level +> telemetry, cloud account, or license-server requirement. Installation may +> contact package repositories to retrieve dependencies. If you purchase through +> Gumroad, Gumroad and its payment providers process checkout, identity, payment +> and tax information under their own policies. +> +> The seller uses order details made available by Gumroad and information you +> voluntarily send for delivery, support, refunds and legally required records. +> Do not send tax IDs, banking information, identity documents, passwords or +> confidential research data to support. We do not need them to troubleshoot +> this local program. Do not post private data in public GitHub issues. +> +> Proposed support-ticket retention is 90 days after closure, except records +> needed for unresolved disputes, security incidents or legal obligations. +> Transaction records are retained as legally required. Contact the seller +> listed on your receipt for privacy requests; identity verification may be +> necessary. No marketing enrollment is required for product delivery. + +Before use, verify actual support tools, providers, storage/access controls, +retention/deletion processes, public contact details, applicable rights and +jurisdictions. Do not claim deletion capabilities or compliance certifications +that have not been implemented. Gumroad's own retention is outside seller control. + +## 4. Validation and provenance statement + +Evidence for commit `58c8ab3dd7aff213e34886d018d828a6dfe0dce6`: + +- [CI run](https://github.com/gmongell/Mathematica_python_TMM_OLEDPixel/actions/runs/38967043233): + eight Windows/Linux Python 3.11-3.14 environments, 56 cases each, + zero failures/errors/skips in parsed JUnit reports. +- Analytic Fresnel interface and Airy film checks, normal-incidence absorbing + characteristic-matrix cross-check, lossless power balance and CLI tests. +- Wheel/source builds, metadata checks, isolated buyer installation, checksums + and runtime dependency audit passed. An advisory scan is not a full security + review and becomes stale as dependencies and advisories change. +- The optional draft-release job was not exercised. No public release is + certified by this document. + +| Evidence or asset | Current status | Remaining action | +|---|---|---| +| Original scripts and contributor rights | Existing repository notices; no complete ownership investigation | Owner reviews history, contributions and any institutional/employer obligations | +| Organic/silver formulas | Inherited placeholders, not material measurements | Trace origin or replace with clearly licensed measured data | +| `tmm`, NumPy, Matplotlib and transitive dependencies | Independently installed; inventory available | Review upstream license/NOTICE obligations for the precise delivery method | +| Scientific reference checks | Bounded plane-wave cases tested | Independent expert review and broader absorbing/oblique benchmarks | +| Images/screenshots | Locally generated preview output available | Select exact-version images with limitations visible | +| Documentation license | Existing CC BY 4.0 notice preserved | Confirm attribution and any adapted third-party content | +| Paid offering | No live listing | Seller verification, approved terms, delivery and checkout tests | + +Existing tests do not establish measured-material accuracy or OLED emission +validity. Do not reuse these results as validation of newly added solver physics. + +## 5. Owner approval record (all pending) + +- [ ] Verify seller legal identity, trade name and authority to sell. +- [ ] Resolve code/contributor/material provenance and license review. +- [ ] Approve limited scientific claims after independent review. +- [ ] Approve price, delivery contents and screenshot set. +- [ ] Approve support capacity, contact, refund and updates policies. +- [ ] Review privacy/legal text against actual operations. +- [ ] Complete seller verification and accept agreements personally. +- [ ] Execute supported test checkout, verify receipt and buyer download. +- [ ] Approve exact release commit, tag/signing process and publication. + +## Official references + +Checked 2026-10-16; platform rules and account-specific requirements can change. + +- [Fees](https://gumroad.com/help/article/66-gumroads-fees) +- [Pricing](https://gumroad.com/pricing) +- [Custom refund policy](https://gumroad.com/help/article/335-custom-refund-policy) +- [Testing a purchase](https://gumroad.com/help/article/62-testing-a-purchase) +- [Terms of service](https://gumroad.com/terms) +- [Privacy policy](https://gumroad.com/privacy) diff --git a/docs/GUMROAD_SELLER_CHECKLIST.md b/docs/GUMROAD_SELLER_CHECKLIST.md new file mode 100644 index 0000000..6f56f37 --- /dev/null +++ b/docs/GUMROAD_SELLER_CHECKLIST.md @@ -0,0 +1,159 @@ +# Gumroad seller verification: guided checklist + +Prepared 2026-10-16. This is operational guidance, not legal/tax advice or +confirmation of approval. Do not put completed identity, tax, address, banking +or document fields into this repository, chat, screenshots or GitHub issues. + +## Current access and next action + +Opening the official [payout settings](https://gumroad.com/settings/payments) +in the available browser redirected to Gumroad login on 2026-10-16. +No authenticated seller profile, verification status or product was accessed. +No form was submitted, account created, agreement accepted or purchase made. + +The next user-only action is to sign in (or register if needed) directly with +Gumroad in a trusted browser, complete any MFA, and review the requested +onboarding steps. Do not send the assistant your password, one-time code, +recovery code, identity document, full account number, SSN or EIN. +Authentication in your usual browser may not authenticate the separate browser +available to this session. + +### Preferred browser: existing Microsoft Edge session + +At the owner's request, the official payout-settings URL was opened using the +running Microsoft Edge installation on 2026-10-16. The tool connection could +not observe/control those existing Edge tabs or confirm their login state. +No Edge process was terminated, profile copied, debugging mode enabled, +extension accessed, or saved password extracted. + +Use Proton Pass normally inside Edge if sign-in is requested. Complete MFA, +identity/tax/bank entry and agreement acceptance yourself on the authenticated +official page. The assistant can interpret non-sensitive instructions/status, +but cannot infer approval or submit invisible forms. + +For help afterward, provide only the non-sensitive field label or error text +and a status such as "identity pending", "bank requested", or "account review". +Omit numbers, personal document images and private dashboard URLs/tokens. + +## 1. Select the truthful account type + +Gumroad's current guidance says: + +- **Individual:** operating under your name, effectively as a sole proprietor, + or without business registration documents. +- **Business:** you have business registration documents; a representative's + personal details and registration documents may also be required. + +The earlier enrollment packets described a sole proprietorship, but this is not +independent evidence of its present legal or registration status. Do not select +"company", claim incorporation, or use a trade name as a substitute for the +verified legal payee unless the actual portal and documents support it. +Use Gumroad support for account-type ambiguity rather than inventing an entity. + +## 2. Prepare privately and enter only in the official portal + +| Portal requirement | Preparation | Who completes it | +|---|---|---| +| Country | Actual supported payout country | Account owner | +| Legal name | Exact name, including middle name where required, matching official records | Account owner | +| Physical address | Valid actual address; Gumroad excludes PO boxes and military addresses in its guidance | Account owner | +| Date of birth/phone | Accurate details; complete any requested verification | Account owner | +| Tax identity | Name/TIN matching IRS records; individual/sole-proprietor treatment depends on actual facts | Account owner, with tax adviser if needed | +| Business evidence | Registration and representative/ownership evidence only if applicable/requested | Authorized representative | +| Photo ID/address evidence | Valid, readable documentation if requested by the secure verification flow | Account owner | +| Payout account | Actual eligible account and payee details matching the selected type | Account owner | +| Agreements | Read and accept personally only with authority | Account owner | + +For US sellers, Gumroad describes EIN use and SSN use for individuals/sole +proprietors without an EIN. Follow the live field instructions and IRS W-9 +instructions for the correct name/TIN combination; do not infer it from a brand. +An EIN does not by itself prove incorporation or remove all personal KYC. + +Gumroad uses Stripe for card processing/bank payouts; the secure flow can request +additional identity verification. PayPal availability depends on the account +and country, not preference alone. Do not switch countries to obtain a payout +method. Gumroad warns country changes can clear saved details and forfeit some +unpaid amounts; contact support first if your country is wrong. + +## 3. Verify each separate status + +- [ ] Login/email and account security complete. +- [ ] Account type and country accurately selected. +- [ ] Legal identity submitted and verification status checked. +- [ ] Tax information submitted and accepted as required. +- [ ] Payout account added and verification status checked. +- [ ] Any business/representative evidence accepted. +- [ ] Any outstanding dashboard warnings resolved. +- [ ] Seller/account review status understood. +- [ ] Product publication eligibility confirmed separately from payouts. + +Do not mark these checks merely because a form was saved. Identity verification, +product eligibility, fraud/account review, and funds actually arriving are +different milestones. Only the platform can confirm acceptance. + +Gumroad's public payout guidance currently describes a standard USD 100 payout +minimum (higher in some locations), holding periods depending on payout schedule, +and possible account review taking 1-3 weeks. These are not guaranteed dates, +account-specific approval criteria, or a reason to create artificial transactions. +Check your actual dashboard and current documentation. + +## 4. Product setup after approval gates + +Keep the product unpublished during preparation. Use +[commercial drafts](COMMERCIAL_DRAFTS.md) as proposed copy, not already approved +terms. Confirm the real product type, supported platforms, price, support +contact, license notices, refund settings and privacy link. + +Before uploading: + +1. Obtain the exact tested wheel/source from the approved build. +2. Verify checksums and audit inventory; include installation/license information. +3. Ensure no customer data, credentials, seller records or development environment + is in the upload. Do not upload `.venv` or a whole workspace backup. +4. Show the alpha and material-data limitations before checkout. +5. Do not promise the not-yet-created OLED solver as a purchase entitlement. +6. Do not enable publication, paid subscriptions or recurring commitments merely + to test a draft. + +## 5. Supported checkout test + +Official Gumroad instructions: remain logged in as the seller, open the product +purchase page and use the built-in **Test card** flow. If the test mode is not +visible, stop and use Gumroad's current supported instructions/support. +The help article also describes a 100% discount-code alternative. + +- Never pay for your own product with your real card to trigger account review + or test payouts. Gumroad explicitly warns this can cause suspension. +- If testing requires publishing or changing product visibility, obtain a + deliberate publication decision first; do not silently make a draft public. +- Check price, tax display, refund-policy visibility, receipt/support link, + filenames, licenses, checksums and download delivery. +- Install the buyer download in a clean environment and run the documented + command. Confirm PNG/CSV/JSON contents against the shipped example. +- A test transaction verifies delivery, not actual settlement or KYC approval. + +## 6. Safe support request template + +> Subject: Seller onboarding/account-type clarification +> +> I am preparing a downloadable scientific/educational Python software product. +> I need help confirming the appropriate account type and the next onboarding +> step for my actual legal status and registration documents. +> The non-sensitive dashboard status is: [insert status only]. +> Please direct me to the secure in-account process for any required evidence. +> I will not send tax identifiers, banking details or identity documents in +> ordinary email. + +Send through Gumroad's official support channel from your authenticated account. +Do not claim production validation or a legal status that has not been verified. + +## Official references + +Checked 2026-10-16; the live portal controls actual requirements. + +- [Payout settings guidance](https://gumroad.com/help/article/260-your-payout-settings-page) +- [Getting paid/account review](https://gumroad.com/help/article/13-getting-paid) +- [Test purchase](https://gumroad.com/help/article/62-testing-a-purchase) +- [Fees](https://gumroad.com/help/article/66-gumroads-fees) +- [Custom refund settings](https://gumroad.com/help/article/335-custom-refund-policy) +- [IRS W-9](https://www.irs.gov/forms-pubs/about-form-w-9) diff --git a/docs/OLED_SOLVER_SPEC.md b/docs/OLED_SOLVER_SPEC.md new file mode 100644 index 0000000..391b14f --- /dev/null +++ b/docs/OLED_SOLVER_SPEC.md @@ -0,0 +1,151 @@ +# OLED dipole-outcoupling solver: first implementation specification + +Prepared 2026-10-16. **Proposed scientific scope, not implemented or validated.** +The plane-wave preview remains separate and must not be relabeled as this solver. +Seller verification is a commercial gate, not a technical prerequisite to +research implementation. A future solver must pass its own scientific gates. + +## Intended first problem + +An electric dipole in a finite, homogeneous, lossless, isotropic emitting layer +of a planar multilayer stack, with passive isotropic surrounding layers. +Use a frequency-domain dyadic Green-function/angular-spectrum formulation and +outgoing-wave boundary conditions. Resolve dipoles parallel and perpendicular +to the interfaces and s/p polarization explicitly. A merely transmitted incident +plane wave is not an emission source and cannot establish outcoupling efficiency. + +Use an explicitly documented time convention, complex-index convention and +square-root branch choice. For example, with `exp(-i omega t)`, passive media +have nonnegative imaginary permittivity and outgoing/decaying longitudinal waves. +Verify that the implementation and any dependency use compatible conventions. + +### Bounded first release + +- Flat laterally infinite layers; no scattering, grating or finite pixel. +- Lossless, positive-index, nonmagnetic emitter host and exterior half-spaces. +- Passive absorbing non-emitting films are allowed with documented + complex permittivity conventions. +- Point dipoles strictly inside the emitter layer, not on an interface. +- No absorbing emitter-host singularity regularization in this first model. +- No exciton transport, electrical injection, thermal effects, stimulated + emission, gain, anisotropy, microcavity lasing or coherence between emitters. +- No claim of full OLED external quantum efficiency from optics alone. + +## Inputs and reproducibility contract + +| Input | Required interpretation | +|---|---| +| Wavelengths | Vacuum wavelength in nm, finite positive ordered samples | +| Stack | Explicit top-to-bottom layer order; first/last are half-spaces | +| Thicknesses | Finite interior thicknesses in nm; positive emitter thickness | +| Material data | Complex n or epsilon convention, units, source, license, interpolation domain | +| Emitter location | Layer index and depth in nm measured from a documented interface | +| Orientation | Vertical fraction in [0,1]; horizontal fraction split equally over in-plane axes | +| Isotropic ensemble | Vertical fraction 1/3, horizontal fraction 2/3 | +| Spectrum | Optional nonnegative intrinsic spectral weights with explicit normalization | +| Intrinsic quantum yield | Optional value in [0,1], with stated reference-medium assumptions | +| Numerical controls | Relative/absolute tolerance, integration strategy and tail-convergence controls | + +Reject nonfinite inputs, invalid orientation/yield, interfaces as dipole +locations, unlicensed/missing material metadata, unsupported host absorption, +and wavelengths outside the measured domain unless explicit, documented +extrapolation is implemented. No silent substitution of placeholder data. +Test numerical interpolation separately from the underlying measurement's validity. + +## Outputs and definitions + +1. Orientation-resolved spectral electromagnetic power relative to a clearly + specified homogeneous-host radiated power for the same dipole. +2. Top and bottom far-field hemispherical powers and angular distributions + with declared angular variables, solid-angle measures and polarization. +3. Total optical absorption per absorbing layer from a physical flux/volume-loss + calculation, not an unexplained residual labeled as absorption. +4. Bound/guided contributions where resolvable, with pole/integration treatment. +5. Numerical balance error, estimated quadrature/tail error and convergence flag. +6. Normalized air-escape fraction and, separately, substrate escape. + In a glass half-space, substrate power is not automatically air outcoupling. + +Let `P_EM` be total source work transferred into electromagnetic channels +(far field, guided and lossy channels) under the specified source convention. +For an accepted passive case, define the optical top-escape fraction as +`P_top_air / P_EM`, not `T_s + T_p`, and not normalized independently per angle. +Do not claim a complete channel decomposition if guided contributions or +high-wavevector lossy tails are unresolved. + +Report raw powers and denominators so users can audit normalization. If a +denominator is zero, integration fails, or errors exceed limits, return an +explicit error/non-converged result; do not clip to [0,1] and report success. + +Only after deriving the rate model may intrinsic yield be combined with the +electromagnetic rate enhancement. For fixed intrinsic nonradiative decay and +`F = P_EM / P_host`, a candidate rate model gives an escape probability per +excitation `q0 * (P_top_air / P_host) / ((1-q0) + q0*F)`. +This requires reference-medium assumptions and special-case tests. It is not +EQE without separately supplied and justified electrical/exciton efficiencies. + +## Numerical architecture + +1. Keep existing preview entry point and outputs backward compatible. +2. Introduce an explicit stack/material/source model, separate from plotting. +3. Implement stable multilayer reflection/response using scattering matrices or + a comparably stable recursion. Do not use exponentially growing transfer + products across thick absorbing/evanescent layers. +4. Derive dipole source terms and Green-function self-interaction for each + orientation. Record equations, references, conventions and units before coding. +5. Integrate in-plane wavevector, separating propagating, evanescent and + guided/pole regions. Split branch points; resolve resonances and near-interface + tails with explicit convergence checks. +6. Use outgoing field flux for far-field channels and physically derived + absorption. Do not count reflected powers as separate emitted populations. +7. Add orientation and emission-depth/spectral ensemble averaging only after + single-dipole normalization is demonstrated. +8. Expose CLI/CSV/JSON only with explicit model version, material provenance, + reference power, tolerances, quadrature diagnostics and limitations. + +## Required validation before commercial claims + +The following are proposed acceptance criteria, not current results: + +| ID | Physical or numerical requirement | Proposed criterion | +|---|---|---| +| O1 | Homogeneous lossless medium, all orientations | Total power / reference = 1 within 1e-6; equal hemispherical total power | +| O2 | Dipole angular patterns | Vertical sin-squared and in-plane azimuthal patterns; integrals within 1e-5 | +| O3 | Symmetric lossless stack, centered emitter | Equal top/bottom power within 1e-6 relative; orientation treated independently | +| O4 | Single-interface dipole reference | Independent analytic/literature comparison across propagating/evanescent regions, agreed tolerance <=1e-4 where reference is well-conditioned | +| O5 | Passive absorbing multilayer | Nonnegative physical powers within error budget; full energy balance residual <=1e-4 of reference power | +| O6 | Integration convergence | Tightened tolerances and enlarged tail domain change integrated outputs <=1e-4 relative with stated absolute floor | +| O7 | Stable thick-film limit | No overflow/NaN; limiting behavior matches independently evaluated boundary cases | +| O8 | Orientation ensemble | Isotropic result equals weighted 1/3 vertical + 2/3 horizontal powers before taking ratios | +| O9 | Near-interface regime | Flag lack of convergence; never report finite-cutoff artifacts as validated efficiency | +| O10 | Material input | Reject out-of-domain wavelengths and unsupported host absorption explicitly | +| O11 | Emission/rate normalization | q0=0 gives zero external photon probability; q0=1 and F=1 reduce to the correct optical limit | +| O12 | Reproducible output | Coordinates, units, layer order, source position, references and error diagnostics round-trip unambiguously | + +Any tolerance adjustment must be justified by numerical/reference error, not +changed merely to obtain passing tests. References and data sets must be +versioned and reviewed for redistribution permissions. + +## Research references and unresolved evidence + +- [S. J. Byrnes, Multilayer optical calculations, arXiv:1603.02720](https://arxiv.org/abs/1603.02720) + supports passive multilayer conventions and plane-wave flux accounting; + it does not by itself supply a complete emitting-dipole model. +- Use the Chance-Prock-Silbey / planar Green-function literature to derive the + source-work and decay-rate expressions; verify the exact source and equations + before implementation. This document does not reproduce a validated derivation. +- Select an independently implemented dipole reference or reproducible published + case with stated stack, constants, orientation, normalization and uncertainties. + The existing `tmm` wrapper cannot serve as an independent dipole reference. +- Obtain measured, licensed optical constants for the intended OLED materials + and measured-stack evidence before any device-specific accuracy claim. + +## Implementation order and sale gates + +**Next technical milestone:** equation/convention derivation plus homogeneous +and single-interface reference kernels/tests (O1, O2, O4), not a full paid solver. +Then stable multilayer source physics, channel balance/convergence, then UI/export. +Keep each stage reviewable and fail closed on unsupported cases. + +Seller onboarding can proceed independently. No publication of an "OLED solver" +offer, pre-order with delivery promises, or reuse of preview test counts as solver +validation until the relevant solver milestones and commercial gates pass. From a40f4070e06d0b3c8f11888eb9d7152ac8f911fa Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 17:52:44 -0400 Subject: [PATCH 4/7] feat: add dipole emission reference kernels and benchmarks Derive and implement normalized homogeneous patterns and single-interface Green-function source-work integrals. Verify image-dipole, dielectric, absorbing and near-field limits with explicit quadrature diagnostics. Preserve the preview CLI and keep finite-stack OLED claims gated. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .github/workflows/plane-wave-preview.yml | 1 + MANIFEST.in | 1 + README.md | 43 ++++ dipole_reference.py | 268 +++++++++++++++++++++ docs/DIPOLE_FORMULATION.md | 153 ++++++++++++ docs/OLED_SOLVER_SPEC.md | 6 + pyproject.toml | 3 +- tests/test_dipole_reference.py | 291 +++++++++++++++++++++++ 8 files changed, 765 insertions(+), 1 deletion(-) create mode 100644 MANIFEST.in create mode 100644 dipole_reference.py create mode 100644 docs/DIPOLE_FORMULATION.md create mode 100644 tests/test_dipole_reference.py diff --git a/.github/workflows/plane-wave-preview.yml b/.github/workflows/plane-wave-preview.yml index bd75a2a..8723f7c 100644 --- a/.github/workflows/plane-wave-preview.yml +++ b/.github/workflows/plane-wave-preview.yml @@ -59,6 +59,7 @@ jobs: "$RUNNER_TEMP/buyer/bin/python" -m pip install dist/*.whl cd "$RUNNER_TEMP" "$RUNNER_TEMP/buyer/bin/mongelli-plane-wave" --wavelength-count 3 --angle-count 3 --output preview.png + "$RUNNER_TEMP/buyer/bin/python" -c "from dipole_reference import SingleInterface; r=SingleInterface(550, 20, 1.0).decay_rates(); assert r.vertical == r.parallel == 1.0" test -s preview.png && test -s preview.csv && test -s preview.json - name: Audit installed runtime dependencies and produce CycloneDX inventory run: | diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..3c2b2c2 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +include docs/DIPOLE_FORMULATION.md diff --git a/README.md b/README.md index 6858c37..4be20a2 100644 --- a/README.md +++ b/README.md @@ -125,6 +125,49 @@ license presentation, final screenshots, seller verification and publication. ## Functions and engineering applications +### Research dipole reference API + +`dipole_reference` provides a separate homogeneous and single-interface +benchmark milestone. It does **not** replace the preview CLI or implement +a finite-stack OLED solver. See [equations, conventions and reference limits](docs/DIPOLE_FORMULATION.md) +and [the larger solver specification](docs/OLED_SOLVER_SPEC.md). + +```python +from dipole_reference import SingleInterface, homogeneous_pattern + +pattern = homogeneous_pattern([0.0, 1.5707963267948966], 0.0, "vertical") +print(pattern.s, pattern.p) # dP / (P_host dOmega), radians, per steradian + +# Synthetic permittivity, not a measured silver data set. +interface = SingleInterface( + wavelength_nm=550.0, distance_nm=20.0, + substrate_epsilon=-5.0 + 0.5j, host_index=1.0, +) +result = interface.decay_rates(atol=1e-8, rtol=1e-7) +print(result.vertical, result.parallel, result.isotropic) +print(result.estimated_error, result.refinement_delta) +``` + +`homogeneous_power` returns macroscopic SI watts for a peak harmonic dipole +amplitude in C m; local-field corrections are not included. +`SingleInterface` requires a lossless real positive host index and positive +dipole distance. Permittivity is relative, with nonnegative imaginary part +under `exp(-i omega t)`; `"pec"` selects the ideal perfect-mirror benchmark. +Lossless nonpositive-permittivity substrates, absorbing hosts, and dimensionless +separations `2*k_host*d < 1e-6` are unsupported and rejected. + +The result is total electromagnetic source work relative to homogeneous-host +power, **not** external efficiency, radiated-only power, or a channel +decomposition. Its propagating/evanescent corrections are integration regions; +in particular host-evanescent fields can transmit into a lossless high-index +substrate. Integration failures raise `DipoleConvergenceError` instead of +returning clipped or success-shaped values. Diagnostic errors are numerical +estimates, not guaranteed bounds or experimental validation. + +The module uses SciPy and is included in the wheel alongside the unchanged +preview API. Tests are in `tests/test_dipole_reference.py`; no measured optical +constants or production OLED claims are supplied. + | Source / interface | Observed function | Inferred application | |---|---|---| | `python_tmm_oled_Green'sfcnApproach.py`: `OLEDOutcouplingTMM` | Layer definitions, aluminum Drude model, interface/propagation matrices, angle/wavelength calculations | Exploring material and thickness sensitivity in OLED stacks [1] | diff --git a/dipole_reference.py b/dipole_reference.py new file mode 100644 index 0000000..bf34ee9 --- /dev/null +++ b/dipole_reference.py @@ -0,0 +1,268 @@ +"""Research dipole reference kernels; not a multilayer OLED outcoupling solver.""" + +from dataclasses import dataclass +from numbers import Real +from typing import Literal + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from scipy.constants import c, epsilon_0 +from scipy.integrate import quad + +Orientation = Literal["vertical", "horizontal_x", "horizontal_y", "horizontal", "isotropic"] +Substrate = complex | Literal["pec"] + + +class DipoleConvergenceError(ArithmeticError): + """The requested source-work accuracy could not be established.""" + + +def _positive(name: str, value: float) -> None: + if not isinstance(value, Real) or not np.isfinite(value) or value <= 0: + raise ValueError(f"{name} must be finite and positive") + + +@dataclass(frozen=True) +class AngularPattern: + """Power per steradian divided by homogeneous-host total power.""" + + s: NDArray[np.float64] + p: NDArray[np.float64] + + @property + def total(self) -> NDArray[np.float64]: + return self.s + self.p + + +def homogeneous_pattern( + theta_rad: ArrayLike, phi_rad: ArrayLike, orientation: Orientation +) -> AngularPattern: + """Broadcast polar/azimuth angles; theta is in [0, pi], phi is periodic.""" + theta, phi = np.broadcast_arrays( + np.asarray(theta_rad, dtype=float), np.asarray(phi_rad, dtype=float) + ) + if not np.all(np.isfinite(theta)) or np.any(theta < 0) or np.any(theta > np.pi): + raise ValueError("Polar angles must be finite and in [0, pi] radians") + if not np.all(np.isfinite(phi)): + raise ValueError("Azimuth angles must be finite radians") + if orientation == "vertical": + s, p = np.zeros_like(theta), np.sin(theta) ** 2 + elif orientation == "horizontal_x": + s, p = np.sin(phi) ** 2, np.cos(theta) ** 2 * np.cos(phi) ** 2 + elif orientation == "horizontal_y": + s, p = np.cos(phi) ** 2, np.cos(theta) ** 2 * np.sin(phi) ** 2 + elif orientation == "horizontal": + s, p = np.full_like(theta, 0.5), 0.5 * np.cos(theta) ** 2 + elif orientation == "isotropic": + s, p = np.full_like(theta, 1 / 3), np.full_like(theta, 1 / 3) + else: + raise ValueError("Unknown dipole orientation") + return AngularPattern(3 / (8 * np.pi) * s, 3 / (8 * np.pi) * p) + + +def homogeneous_power(wavelength_nm: float, host_index: float, dipole_c_m: complex) -> float: + """Macroscopic SI watts for peak harmonic dipole amplitude; no local-field factor.""" + _positive("Vacuum wavelength (nm)", wavelength_nm) + _positive("Host refractive index", host_index) + if not np.isfinite(dipole_c_m): + raise ValueError("Dipole amplitude must be finite (C m)") + with np.errstate(over="raise", divide="raise", invalid="raise"): + try: + omega = 2 * np.pi * c / (np.float64(wavelength_nm) * 1e-9) + power = host_index * omega**4 * np.abs(np.complex128(dipole_c_m)) ** 2 / (12 * np.pi * epsilon_0 * c**3) + except FloatingPointError as error: + raise ValueError("Homogeneous power overflowed; check input units") from error + if not np.isfinite(power): + raise ValueError("Homogeneous power overflowed; check input units") + return float(power) + + +@dataclass(frozen=True) +class DecayResult: + """Total electromagnetic source-work ratios, not escape efficiencies. + + Tuple diagnostics are ordered (vertical, parallel). Corrections can be + negative and are integration regions, not individually positive channels. + """ + + vertical: float + parallel: float + propagating_correction: tuple[float, float] + evanescent_correction: tuple[float, float] + estimated_error: tuple[float, float] + refinement_delta: tuple[float, float] + tail_magnitude: tuple[float, float] + evaluations: int + + def oriented(self, vertical_fraction: float) -> float: + if not np.isfinite(vertical_fraction) or not 0 <= vertical_fraction <= 1: + raise ValueError("Vertical dipole fraction must be finite and in [0, 1]") + return vertical_fraction * self.vertical + (1 - vertical_fraction) * self.parallel + + @property + def isotropic(self) -> float: + return self.oriented(1 / 3) + + +@dataclass(frozen=True) +class SingleInterface: + """Dipole at z=d>0 above a substrate at z<0; host is lossless.""" + + wavelength_nm: float + distance_nm: float + substrate_epsilon: Substrate + host_index: float = 1.0 + + def __post_init__(self) -> None: + _positive("Vacuum wavelength (nm)", self.wavelength_nm) + _positive("Dipole-interface distance (nm)", self.distance_nm) + _positive("Lossless host refractive index", self.host_index) + if not np.isfinite(self.phase_distance) or self.phase_distance < 1e-6: + raise ValueError("Dimensionless separation 2*k_host*d must be finite and >= 1e-6") + if isinstance(self.substrate_epsilon, str): + if self.substrate_epsilon != "pec": + raise ValueError("The only ideal-boundary sentinel is 'pec'") + else: + epsilon = complex(self.substrate_epsilon) + if not np.isfinite(epsilon) or epsilon.imag < 0: + raise ValueError("Substrate permittivity must be finite and passive (Im >= 0)") + if epsilon.imag == 0 and epsilon.real <= 0: + raise ValueError("Lossless nonpositive-epsilon substrates require unsupported pole/ENZ treatment") + + @property + def phase_distance(self) -> float: + return 4 * np.pi * self.host_index * self.distance_nm / self.wavelength_nm + + def _reflection(self, w: complex) -> tuple[complex, complex]: + if self.substrate_epsilon == "pec": + return -1 + 0j, 1 + 0j + relative = complex(self.substrate_epsilon) / self.host_index**2 + if relative == 1: + return 0j, 0j + v = np.sqrt(relative - 1 + w * w + 0j) + return (w - v) / (w + v), (relative * w - v) / (relative * w + v) + + def spectral_correction(self, q: float) -> tuple[float, float]: + """Vertical/parallel d(F-1)/dq; rejects the singular host light line q=1.""" + if not np.isfinite(q) or q < 0 or q == 1: + raise ValueError("q=k_parallel/k_host must be finite, nonnegative, and not 1") + w = np.sqrt(1 - q * q + 0j) + rs, rp = self._reflection(w) + phase = np.exp(1j * self.phase_distance * w) + values = ( + 1.5 * np.real(q**3 / w * rp * phase), + 0.75 * np.real(q / w * (rs - w * w * rp) * phase), + ) + if not np.all(np.isfinite(values)): + raise DipoleConvergenceError("Nonfinite raw spectral kernel; use a supported q range") + return float(values[0]), float(values[1]) + + def _breakpoints(self) -> tuple[list[float], list[float]]: + prop: list[float] = [] + evan: list[float] = [] + if self.substrate_epsilon == "pec": + return prop, evan + relative = complex(self.substrate_epsilon) / self.host_index**2 + if relative.imag == 0: + if 0 < relative.real < 1: + prop.append(float(np.sqrt(1 - relative.real))) + elif relative.real > 1: + evan.append(float(self.phase_distance * np.sqrt(relative.real - 1))) + elif relative.real < -1: + pole = self.phase_distance * np.sqrt(-1 / (relative + 1)) + center, width = float(pole.real), abs(float(pole.imag)) + offsets = (-32768, -4096, -512, -64, -8, -1, 0, 1, 8, 64, 512, 4096, 32768) + evan.extend(center + offset * width for offset in offsets) + elif relative.imag > 0: + # Resolve a weakly lossy dielectric's smoothed substrate light line. + if relative.real < 1: + branch = np.sqrt(1 - relative) + prop.extend( + point for offset in (-8, -1, 0, 1, 8) + if 0 < (point := float(branch.real + offset * abs(branch.imag))) < 1 + ) + elif relative.real > 1: + branch = self.phase_distance * np.sqrt(relative - 1) + evan.extend(float(branch.real + offset * abs(branch.imag)) for offset in (-8, -1, 0, 1, 8)) + return sorted(set(prop)), sorted(set(point for point in evan if point > 0)) + + def decay_rates( + self, *, atol: float = 1e-8, rtol: float = 1e-7, max_subdivisions: int = 300 + ) -> DecayResult: + """Integrate source work, repeating transformed infinite-tail quadrature. + + Raises DipoleConvergenceError on failed error estimates/refinement. + No quadrature estimate is a rigorous bound on all possible resonances. + """ + if not 1e-12 <= atol <= 1e-2 or not 1e-12 <= rtol <= 1e-2: + raise ValueError("atol and rtol must be finite and in [1e-12, 1e-2]") + if isinstance(max_subdivisions, bool) or not isinstance(max_subdivisions, int) or not 1 <= max_subdivisions <= 10000: + raise ValueError("max_subdivisions must be an integer in [1, 10000]") + prop_points, evan_points = self._breakpoints() + a = self.phase_distance + evaluations = 0 + + def integrate(function, boundaries, tolerance): + nonlocal evaluations + total = error = 0.0 + budget = atol * tolerance / (len(boundaries) - 1) + for left, right in zip(boundaries[:-1], boundaries[1:]): + outcome = quad( + function, left, right, epsabs=budget, + epsrel=rtol * tolerance, limit=max_subdivisions, full_output=1, + ) + value, estimate, info = outcome[:3] + evaluations += info["neval"] + if len(outcome) != 3 or not np.isfinite(value) or not np.isfinite(estimate): + message = outcome[3] if len(outcome) > 3 else "nonfinite quadrature" + raise DipoleConvergenceError(f"Single-interface integration failed: {message}") + total += value + error += estimate + return total, error + + corrections, errors, deltas, tails = [], [], [], [] + for vertical in (True, False): + factor = 1.5 if vertical else 0.75 + + def prop(u): + rs, rp = self._reflection(complex(u)) + amplitude = (1 - u * u) * rp if vertical else rs - u * u * rp + return float(factor * np.real(amplitude * np.exp(1j * a * u))) + + def evan(x): + if x > 700: + return 0.0 # Beyond floating-point exponential support at the requested accuracy. + t = x / a + rs, rp = self._reflection(1j * t) + amplitude = (1 + t * t) * rp if vertical else rs + t * t * rp + return float(factor / a * np.exp(-x) * np.imag(amplitude)) + + passes = [] + for split, tolerance in ((32.0, 0.05), (64.0, 0.0125)): + p, ep = integrate(prop, [0.0, *prop_points, 1.0], tolerance) + points = sorted(set([0.0, split, *evan_points])) + e, ee = integrate(evan, [*points, np.inf], tolerance) + tail, _ = integrate(evan, [points[-1], np.inf], tolerance) + passes.append((p, e, ep + ee, abs(tail))) + p, e, error, tail = passes[-1] + rate = 1 + p + e + delta = abs(p - passes[0][0]) + abs(e - passes[0][1]) + budget = max(atol, rtol * abs(rate)) + if not np.isfinite(rate) or rate < 0 or max(error, delta) > budget: + raise DipoleConvergenceError( + f"Source-work accuracy not established: rate={rate}, error={error}, refinement={delta}" + ) + corrections.append((p, e)) + errors.append(error) + deltas.append(delta) + tails.append(tail) + return DecayResult( + vertical=1 + sum(corrections[0]), + parallel=1 + sum(corrections[1]), + propagating_correction=(corrections[0][0], corrections[1][0]), + evanescent_correction=(corrections[0][1], corrections[1][1]), + estimated_error=(errors[0], errors[1]), + refinement_delta=(deltas[0], deltas[1]), + tail_magnitude=(tails[0], tails[1]), + evaluations=evaluations, + ) diff --git a/docs/DIPOLE_FORMULATION.md b/docs/DIPOLE_FORMULATION.md new file mode 100644 index 0000000..8f4a5a5 --- /dev/null +++ b/docs/DIPOLE_FORMULATION.md @@ -0,0 +1,153 @@ +# Dipole reference kernels: conventions and derivation + +Research milestone, 2026-10-16. Implemented in `dipole_reference.py`. +These are homogeneous and **one-interface** reference kernels, not a finite +emitting-layer/multilayer OLED solver, material database, or sale-ready product. + +## Geometry and normalization + +The interface is z=0. A prescribed harmonic electric dipole is at z=d>0 in +a lossless isotropic nonmagnetic host with real n1>0, epsilon1=n1^2. +The substrate is the z<0 half-space, with relative permittivity epsilon2. +Use exp(-i omega t), SI dipole amplitude p (C m), vacuum wavelength lambda0 +(nm at the API), and k1=2 pi n1/lambda0 (inverse nm for phase calculations). +The substrate is passive: Im(epsilon2)>=0. The API also accepts the explicit +ideal-boundary sentinel `"pec"` for a perfect electric conductor. + +Homogeneous macroscopic radiated power, for peak (not RMS) dipole amplitude: + + P_host = n1 * omega^4 * |p|^2 / (12 pi epsilon0 c^3) + +No microscopic local-field correction is included. The normalized angular +density is dP/(P_host dOmega), with dOmega=sin(theta) dtheta dphi. +For a unit dipole axis e and observation direction r_hat: + + density = 3/(8 pi) * (1 - |e . r_hat|^2) + +Project e onto e_phi (s) and e_theta (p) before squaring. In particular: + +| Dipole | s density / [3/(8 pi)] | p density / [3/(8 pi)] | +|---|---|---| +| vertical z | 0 | sin(theta)^2 | +| horizontal x | sin(phi)^2 | cos(theta)^2 cos(phi)^2 | +| horizontal y | cos(phi)^2 | cos(theta)^2 sin(phi)^2 | +| random in-plane | 1/2 | cos(theta)^2/2 | +| isotropic ensemble | 1/3 | 1/3 | + +Each total density integrates to 1 over the sphere and 1/2 in each hemisphere. +An isotropic ensemble averages **powers**, not amplitudes: vertical weight +1/3 plus random-in-plane weight 2/3. + +## Reflected Green tensor and source work + +Let q=k_parallel/k1, w=sqrt(1-q^2), and +v=sqrt(epsilon2/epsilon1-q^2). Select the outgoing/decaying branch: +positive real root for propagation, positive imaginary root for evanescence. +The p coefficient here uses the standard electric-field Fresnel convention: + + r_s = (w-v)/(w+v) + r_p = (epsilon2*w-epsilon1*v)/(epsilon2*w+epsilon1*v) + +For `"pec"`, r_s=-1 and r_p=+1. Identical media give both coefficients zero. +The total electromagnetic source-work enhancement relative to P_host is: + + F_vertical = 1 + (3/2) Re integral_0^infinity + [q^3/w * r_p * exp(2 i k1 d w)] dq + + F_parallel = 1 + (3/4) Re integral_0^infinity + [q/w * (r_s - w^2 r_p) * exp(2 i k1 d w)] dq + +The **minus** before w^2*r_p in the parallel expression is essential. It +comes from the reflected p dyad, not from a polarization-intensity average. +These expressions follow the reflected dyadic Green-function contribution +to source work (Im[p* . E_reflected]) with homogeneous normalization. +The in-plane orientations are equivalent after azimuthal integration. + +These rates exclude intrinsic emitter nonradiative decay. For an absorbing +substrate they include work dissipated in that substrate. They are NOT +far-field escape fractions, photon yield, or external quantum efficiency. +In particular, a host-evanescent component can propagate in a higher-index +lossless substrate; "evanescent correction" does not mean "absorbed power". + +## Integration without the host light-line singularity + +Write a=2 k1 d. In the q<1 region substitute u=sqrt(1-q^2): + + Delta_vertical_prop = (3/2) integral_0^1 + (1-u^2) Re[r_p exp(i a u)] du + Delta_parallel_prop = (3/4) integral_0^1 + Re[(r_s-u^2*r_p) exp(i a u)] du + +In q>1 substitute t=sqrt(q^2-1), then x=a*t: + + Delta_vertical_evan = (3/2)/a integral_0^infinity + (1+(x/a)^2) exp(-x) Im[r_p] dx + Delta_parallel_evan = (3/4)/a integral_0^infinity + exp(-x) Im[r_s+(x/a)^2*r_p] dx + +This removes the integrable q=1 singularity and scales the decaying tail. +For a lossless dielectric, split the substrate light line explicitly: +u=sqrt(1-epsilon2/epsilon1) or x=a*sqrt(epsilon2/epsilon1-1). +For a lossy metal, also split near the complex surface-mode location +t=sqrt(-1/(epsilon2/epsilon1+1)), including its width. Strictly lossless +nonpositive-epsilon substrates are rejected: real-axis poles and ENZ limits +need an explicit limiting-absorption/pole prescription not supplied here. + +Weakly absorbing dielectric substrate light lines are also split near their +complex branch location. Surface-mode splits span geometrically increasing +width offsets so narrow resonances and their wings are not left in one large +integration interval. Extremely narrow poles can still fail the accuracy gate. + +The full infinite evanescent tail is integrated, not truncated at a chosen +wavevector. Split it initially at x=32, then repeat with x=64 and four-times +tighter quadrature tolerances. Return estimated quadrature error, difference +between refinements, the absolute signed tail contribution, and evaluation count. +The tail diagnostic is not a rigorous error bound. Reject QUADPACK failures, +nonfinite results, negative total rates, or errors/refinement changes above +max(atol, rtol*|F|); never clip a failed rate into a plausible answer. +Very small dimensionless separation a<1e-6 is explicitly unsupported. +Even above that numerical boundary, continuum optical constants may be invalid +at molecular distances. Quadrature estimates are not a proof of convergence +for arbitrarily sharp resonances; independent benchmarks remain necessary. + +The raw q-space kernel intentionally rejects q=1 (the singular point). +Use the transformed integrated API, not sampling at the light line. + +## Analytic benchmark: ideal mirror + +Method of images gives, with a=2 k1 d: + + F_vertical = 1 + 3*(sin(a)-a*cos(a))/a^3 + F_parallel = 1 - (3/2)*(sin(a)/a + cos(a)/a^2 - sin(a)/a^3) + +At small a, the numerically stable limits are +F_vertical=2-a^2/10+a^4/280+... and +F_parallel=a^2/5-3*a^4/280+.... +At large separation both tend to 1. These independently derived image-dipole +results detect wrong signs, missing prefactors and phase/unit mistakes. + +For a lossy half-space, the quasistatic leading evanescent corrections are: + + Delta_vertical ~ 3/a^3 * Im[(epsilon2-epsilon1)/(epsilon2+epsilon1)] + Delta_parallel ~ (3/2)/a^3 * Im[(epsilon2-epsilon1)/(epsilon2+epsilon1)] + +This is a mathematical asymptote, not a material/device-validity claim. + +## Sources and validation limits + +- [Auguie, planar: decay-rate theory](https://nano-optics.ac.nz/planar/articles/decay_rates_theory.html), + accessed 2026-10-16. Explicit orientation formulas and variable substitutions; + the implementation here is newly written from the equations, not copied code. +- [Novotny and Hecht, Principles of Nano-Optics, chapter 10](https://www.cambridge.org/core/books/principles-of-nanooptics/dipole-emission-near-planar-interfaces/059D85E0BA32D70B1B9FCB83A4ABB62A). + Broader planar-interface Green-function framework. This citation is not an + assertion that all chapter cases are implemented or independently validated. +- [Byrnes, Multilayer optical calculations](https://arxiv.org/abs/1603.02720). + Background on passive-media conventions and Fresnel power accounting. + +Tests compare mirror image-dipole limits, full-sphere homogeneous normalization, +an alternative theta/hyperbolic integration for synthetic dielectric/absorbing +interfaces, dielectric branch cases and a quasistatic asymptote. These are +bounded reference checks, not independent experimental validation. No measured +material constants or redistributed external reference data are supplied. +Finite stacks, emitter-layer feedback, layer-resolved absorption, guided-mode +accounting and air/substrate escape still require subsequent implementation. diff --git a/docs/OLED_SOLVER_SPEC.md b/docs/OLED_SOLVER_SPEC.md index 391b14f..85c895b 100644 --- a/docs/OLED_SOLVER_SPEC.md +++ b/docs/OLED_SOLVER_SPEC.md @@ -5,6 +5,12 @@ The plane-wave preview remains separate and must not be relabeled as this solver Seller verification is a commercial gate, not a technical prerequisite to research implementation. A future solver must pass its own scientific gates. +**First reference milestone implemented:** `dipole_reference.py` contains +homogeneous angular/power kernels and single-interface total source-work +integration, with [documented equations and benchmark scope](DIPOLE_FORMULATION.md). +This partially addresses O1/O2/O4/O6/O8/O9 only for the stated reference geometry; +it does not satisfy the complete finite-multilayer acceptance table below. + ## Intended first problem An electric dipole in a finite, homogeneous, lossless, isotropic emitting layer diff --git a/pyproject.toml b/pyproject.toml index 3794ce0..87e3695 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,6 +15,7 @@ dependencies = [ "numpy>=1.24,<3", "matplotlib>=3.7,<4", "tmm>=0.2.0,<0.3", + "scipy>=1.11,<2", ] classifiers = [ "Development Status :: 3 - Alpha", @@ -32,7 +33,7 @@ mongelli-plane-wave = "tmm_OLED9:main" Repository = "https://github.com/gmongell/Mathematica_python_TMM_OLEDPixel" [tool.setuptools] -py-modules = ["tmm_OLED9"] +py-modules = ["tmm_OLED9", "dipole_reference"] [tool.pytest.ini_options] testpaths = ["tests"] diff --git a/tests/test_dipole_reference.py b/tests/test_dipole_reference.py new file mode 100644 index 0000000..bdbc27e --- /dev/null +++ b/tests/test_dipole_reference.py @@ -0,0 +1,291 @@ +import numpy as np +import pytest +from scipy.constants import c, epsilon_0 +from scipy.integrate import quad + +from dipole_reference import ( + DipoleConvergenceError, + SingleInterface, + homogeneous_pattern, + homogeneous_power, +) + + +@pytest.mark.parametrize("orientation", [ + "vertical", "horizontal_x", "horizontal_y", "horizontal", "isotropic", +]) +def test_homogeneous_sphere_and_hemisphere_normalization(orientation): + nodes, weights = np.polynomial.legendre.leggauss(40) + phi = np.linspace(0, 2 * np.pi, 80, endpoint=False) + for lo, hi, expected in [(-1, 1, 1), (0, 1, 0.5), (-1, 0, 0.5)]: + cosine = lo + (nodes + 1) * (hi - lo) / 2 + density = homogeneous_pattern(np.arccos(cosine[:, None]), phi, orientation) + integral = np.sum(weights[:, None] * density.total) * (hi - lo) / 2 * 2 * np.pi / len(phi) + assert integral == pytest.approx(expected, abs=1e-12) + assert np.all(density.s >= 0) and np.all(density.p >= 0) + + +def test_homogeneous_directional_patterns_and_polarization(): + normalization = 3 / (8 * np.pi) + vertical = homogeneous_pattern([0, np.pi / 2, np.pi], 0, "vertical") + np.testing.assert_allclose(vertical.total, [0, normalization, 0], atol=1e-15) + np.testing.assert_array_equal(vertical.s, [0, 0, 0]) + x = homogeneous_pattern(np.pi / 2, [0, np.pi / 2], "horizontal_x") + y = homogeneous_pattern(np.pi / 2, [0, np.pi / 2], "horizontal_y") + np.testing.assert_allclose(x.total, [0, normalization], atol=1e-15) + np.testing.assert_allclose(y.total, [normalization, 0], atol=1e-15) + np.testing.assert_allclose(x.p, [0, 0], atol=1e-15) + isotropic = homogeneous_pattern([0, 0.7, np.pi], [1, 2, 3], "isotropic") + np.testing.assert_allclose(isotropic.s, 1 / (8 * np.pi), atol=1e-15) + np.testing.assert_allclose(isotropic.p, 1 / (8 * np.pi), atol=1e-15) + + +def test_isotropic_pattern_is_orientation_power_average(): + theta, phi = np.linspace(0, np.pi, 17)[:, None], np.linspace(0, 2 * np.pi, 19) + vertical = homogeneous_pattern(theta, phi, "vertical") + horizontal = homogeneous_pattern(theta, phi, "horizontal") + x = homogeneous_pattern(theta, phi, "horizontal_x") + y = homogeneous_pattern(theta, phi, "horizontal_y") + isotropic = homogeneous_pattern(theta, phi, "isotropic") + for part in ("s", "p"): + np.testing.assert_allclose(getattr(horizontal, part), (getattr(x, part) + getattr(y, part)) / 2, atol=1e-15) + np.testing.assert_allclose( + getattr(isotropic, part), (getattr(vertical, part) + 2 * getattr(horizontal, part)) / 3, + atol=1e-15, + ) + + +def test_homogeneous_power_si_normalization_and_scaling(): + wavelength_m, amplitude = 550e-9, 1e-29 + # Equivalent k0^4*c expression checks wavelength conversion and harmonic amplitude. + expected = (2 * np.pi / wavelength_m) ** 4 * c * amplitude**2 / (12 * np.pi * epsilon_0) + vacuum = homogeneous_power(550, 1, amplitude) + assert vacuum == pytest.approx(expected, rel=1e-13) + assert homogeneous_power(550, 1.5, amplitude) == pytest.approx(1.5 * vacuum) + assert homogeneous_power(1100, 1, 2j * amplitude) == pytest.approx(vacuum / 4) + assert homogeneous_power(550, 1, 0) == 0 + + +@pytest.mark.parametrize("host_index", [1.0, 1.7]) +def test_matched_interface_has_unit_rate_and_zero_spectrum(host_index): + interface = SingleInterface(550, 37, host_index**2, host_index) + for q in [0, 0.5, 0.999, 1.001, 2, 10]: + assert interface.spectral_correction(q) == (0, 0) + result = interface.decay_rates() + assert result.vertical == result.parallel == result.isotropic == 1 + assert result.propagating_correction == result.evanescent_correction == (0, 0) + assert result.estimated_error == result.tail_magnitude == (0, 0) + + +@pytest.mark.parametrize("a", [0.001, 0.1, 1.0, 4.0, 20.0, 100.0]) +def test_perfect_mirror_rates_match_image_dipole(a): + interface = SingleInterface(600, a * 600 / (4 * np.pi * 1.5), "pec", 1.5) + result = interface.decay_rates(atol=1e-10, rtol=1e-9) + if a < 0.01: + expected_vertical = 2 - a**2 / 10 + a**4 / 280 + expected_parallel = a**2 / 5 - 3 * a**4 / 280 + else: + expected_vertical = 1 + 3 * (np.sin(a) - a * np.cos(a)) / a**3 + expected_parallel = 1 - 1.5 * (np.sin(a) / a + np.cos(a) / a**2 - np.sin(a) / a**3) + assert result.vertical == pytest.approx(expected_vertical, abs=1e-10) + assert result.parallel == pytest.approx(expected_parallel, abs=1e-10) + assert result.evanescent_correction == (0, 0) + assert result.evaluations > 0 + + +def test_mirror_raw_spectrum_has_correct_parallel_sign(): + interface = SingleInterface(600, 50, "pec") + q = 0.6 + w = np.sqrt(1 - q**2) + phase_real = np.cos(4 * np.pi * 50 / 600 * w) + vertical, parallel = interface.spectral_correction(q) + assert vertical == pytest.approx(1.5 * q**3 / w * phase_real) + assert parallel == pytest.approx(-0.75 * q / w * (1 + w**2) * phase_real) + + +def _alternative_angular_integral(epsilon, distance, wavelength, host): + """Different theta/cosh coordinates and absolute wavevectors for reference.""" + k0, k1 = 2 * np.pi / wavelength, 2 * np.pi * host / wavelength + + def reflection(q): + kz1 = k1 * np.sqrt(1 - q * q + 0j) + kz2 = np.sqrt(k0 * k0 * epsilon - (q * k1) ** 2 + 0j) + rs = (kz1 - kz2) / (kz1 + kz2) + rp = (epsilon * kz1 - host**2 * kz2) / (epsilon * kz1 + host**2 * kz2) + return rs, rp + + values = [] + for vertical in [True, False]: + def angular(theta): + rs, rp = reflection(np.sin(theta)) + cosine = np.cos(theta) + amplitude = np.sin(theta) ** 3 * rp if vertical else np.sin(theta) * (rs - cosine**2 * rp) + return np.real(amplitude * np.exp(2j * k1 * distance * cosine)) + + def hyperbolic(t): + rs, rp = reflection(np.cosh(t)) + sinh = np.sinh(t) + amplitude = np.cosh(t) ** 3 * rp if vertical else np.cosh(t) * (rs + sinh**2 * rp) + return np.imag(amplitude) * np.exp(-2 * k1 * distance * sinh) + + prop_points, evan_points = [], [] + if complex(epsilon).imag == 0: + ratio = float(np.real(epsilon)) / host**2 + if ratio < 1: + prop_points = [np.arcsin(np.sqrt(ratio))] + elif ratio > 1: + evan_points = [np.arccosh(np.sqrt(ratio))] + prop = quad(angular, 0, np.pi / 2, points=prop_points, epsabs=1e-11, epsrel=1e-11, limit=300)[0] + upper = np.arcsinh(80 / (2 * k1 * distance)) + evan = quad(hyperbolic, 0, upper, points=[p for p in evan_points if p < upper], + epsabs=1e-11, epsrel=1e-11, limit=300)[0] + values.append(1 + (1.5 if vertical else 0.75) * (prop + evan)) + return values + + +@pytest.mark.parametrize("epsilon,host", [ + (2.25, 1.0), (1.0, 1.5), (2.25 + 0.4j, 1.0), (-5 + 0.5j, 1.0), +]) +def test_interface_matches_alternative_angular_integration(epsilon, host): + result = SingleInterface(550, 45, epsilon, host).decay_rates(atol=1e-10, rtol=1e-9) + expected = _alternative_angular_integral(epsilon, 45, 550, host) + np.testing.assert_allclose([result.vertical, result.parallel], expected, rtol=1e-8, atol=1e-9) + + +def test_high_index_lossless_substrate_has_host_evanescent_source_work(): + result = SingleInterface(550, 40, 2.25).decay_rates() + assert result.evanescent_correction[0] > 0.01 + assert result.evanescent_correction[1] > 0.01 + # This contribution is transmission into a dielectric, not substrate absorption. + assert result.vertical > 0 and result.parallel > 0 + + +def test_lossy_quasistatic_limit_and_orientation_ratio(): + epsilon = -5 + 0.5j + interface = SingleInterface(550, 0.05, epsilon) + result = interface.decay_rates(atol=1e-8, rtol=1e-7) + a = 4 * np.pi * 0.05 / 550 + reflection_imag = np.imag((epsilon - 1) / (epsilon + 1)) + assert (result.vertical - 1) * a**3 / 3 == pytest.approx(reflection_imag, rel=2e-5) + assert (result.parallel - 1) * a**3 / 1.5 == pytest.approx(reflection_imag, rel=2e-5) + + +def test_refinement_and_diagnostics_for_lossy_interface(): + interface = SingleInterface(550, 12, -3 + 0.2j) + loose = interface.decay_rates() + tight = interface.decay_rates(atol=1e-10, rtol=1e-9) + np.testing.assert_allclose([loose.vertical, loose.parallel], [tight.vertical, tight.parallel], rtol=1e-7) + for i, rate in enumerate([tight.vertical, tight.parallel]): + assert tight.estimated_error[i] <= max(1e-10, 1e-9 * rate) + assert tight.refinement_delta[i] <= max(1e-10, 1e-9 * rate) + assert 0 <= tight.tail_magnitude[i] < 1e-12 + assert rate == pytest.approx(1 + tight.propagating_correction[i] + tight.evanescent_correction[i]) + assert tight.oriented(0) == tight.parallel + assert tight.oriented(1) == tight.vertical + assert tight.isotropic == pytest.approx((tight.vertical + 2 * tight.parallel) / 3) + + +def test_dimensionless_interface_scaling_preserves_source_work(): + baseline = SingleInterface(550, 20, -5 + 0.5j).decay_rates() + scaled_lengths = SingleInterface(1100, 40, -5 + 0.5j).decay_rates() + scaled_medium = SingleInterface(1100, 20, 4 * (-5 + 0.5j), 2).decay_rates() + for other in [scaled_lengths, scaled_medium]: + np.testing.assert_allclose( + [other.vertical, other.parallel], [baseline.vertical, baseline.parallel], rtol=1e-12, + ) + + +@pytest.mark.parametrize("epsilon,host", [(2.25, 1.0), (1.0, 1.5)]) +def test_weak_dielectric_loss_approaches_lossless_branch_limit(epsilon, host): + lossless = SingleInterface(550, 20, epsilon, host).decay_rates() + weak_loss = SingleInterface(550, 20, epsilon + 1e-7j, host).decay_rates() + np.testing.assert_allclose( + [weak_loss.vertical, weak_loss.parallel], [lossless.vertical, lossless.parallel], + rtol=1e-5, atol=1e-7, + ) + + +def test_weak_loss_surface_resonance_is_resolved(): + interface = SingleInterface(550, 20, -5 + 1e-5j) + result = interface.decay_rates() + refined = interface.decay_rates(atol=1e-10, rtol=1e-9) + np.testing.assert_allclose( + [result.vertical, result.parallel], [refined.vertical, refined.parallel], rtol=1e-7, + ) + assert result.evanescent_correction[0] > 0.1 + assert result.evanescent_correction[1] > 0.01 + + +def test_far_mirror_tends_to_homogeneous_rate(): + result = SingleInterface(550, 100 * 550 / (4 * np.pi), "pec").decay_rates() + assert result.vertical == pytest.approx(1, abs=4e-4) + assert result.parallel == pytest.approx(1, abs=0.016) + + +def test_lossy_raw_evanescent_spectrum_is_nonzero_and_positive(): + interface = SingleInterface(550, 20, -5 + 0.5j) + vertical, parallel = interface.spectral_correction(2) + assert vertical > 0.01 + assert parallel > 0.01 + # At high q the p-dipole near-field terms approach a 2:1 ratio. + vertical, parallel = interface.spectral_correction(20) + assert vertical / parallel == pytest.approx(2, rel=0.01) + + +def test_underresolved_quadrature_fails_explicitly(): + interface = SingleInterface(550, 2000, "pec") + with pytest.raises(DipoleConvergenceError, match="integration failed"): + interface.decay_rates(max_subdivisions=1) + + +@pytest.mark.parametrize("wavelength,distance,epsilon,host", [ + (0, 20, 2.25, 1), (np.nan, 20, 2.25, 1), (np.inf, 20, 2.25, 1), + (550, 0, 2.25, 1), (550, -1, 2.25, 1), (550, np.nan, 2.25, 1), + (550, 1e-12, 2.25, 1), (550, 20, 2.25, 0), (550, 20, 2.25, np.inf), + (550, 20, -1, 1), (550, 20, 0, 1), (550, 20, 2.25 - 0.1j, 1), + (550, 20, complex(np.nan, 1), 1), (550, 20, complex(1, np.inf), 1), + (550, 20, "mirror", 1), + (550, 20, 2.25, 1 + 0.1j), +]) +def test_invalid_interface_inputs(wavelength, distance, epsilon, host): + with pytest.raises(ValueError): + SingleInterface(wavelength, distance, epsilon, host) + + +@pytest.mark.parametrize("q", [-1, 1, np.nan, np.inf]) +def test_invalid_spectral_coordinates(q): + with pytest.raises(ValueError, match="q="): + SingleInterface(550, 20, 2.25).spectral_correction(q) + + +@pytest.mark.parametrize("kwargs", [ + {"atol": 0}, {"rtol": np.nan}, {"rtol": np.inf}, {"atol": 1e-15}, + {"max_subdivisions": 0}, {"max_subdivisions": 2.5}, {"max_subdivisions": True}, +]) +def test_invalid_integration_controls(kwargs): + with pytest.raises(ValueError): + SingleInterface(550, 20, 2.25).decay_rates(**kwargs) + + +@pytest.mark.parametrize("theta,phi,orientation", [ + (-1, 0, "vertical"), (np.pi + 0.01, 0, "vertical"), (np.nan, 0, "vertical"), + (0, np.inf, "vertical"), (0, 0, "unknown"), +]) +def test_invalid_pattern_inputs(theta, phi, orientation): + with pytest.raises(ValueError): + homogeneous_pattern(theta, phi, orientation) + + +@pytest.mark.parametrize("wavelength,host,amplitude", [ + (0, 1, 1e-29), (550, 0, 1e-29), (550, 1, np.nan), (550, 1, np.inf), + (550, 1 + 0.1j, 1e-29), (1e-100, 1, 1e100), +]) +def test_invalid_homogeneous_power_inputs(wavelength, host, amplitude): + with pytest.raises(ValueError): + homogeneous_power(wavelength, host, amplitude) + + +@pytest.mark.parametrize("fraction", [-0.01, 1.01, np.nan, np.inf]) +def test_invalid_orientation_fraction(fraction): + with pytest.raises(ValueError): + SingleInterface(550, 20, "pec").decay_rates().oriented(fraction) From 939b689baccb3860d519c5d6743f12258f4a2423 Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 18:22:33 -0400 Subject: [PATCH 5/7] feat: add spectral data, EQE composition and ensemble planning Track native dispersive material provenance and compose explicit optical-rate EQE. Compare incoherent orientations and estimate planar versus finite-pixel resources without claiming unimplemented source, GPU or FDTD physics. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .github/workflows/plane-wave-preview.yml | 1 + MANIFEST.in | 2 + README.md | 15 + docs/ENSEMBLE_RESOURCE_ESTIMATES.md | 169 +++++++++ docs/SPECTRAL_EQE_AND_DATA.md | 210 ++++++++++ pyproject.toml | 4 +- simulation_resources.py | 266 +++++++++++++ spectral_eqe.py | 463 +++++++++++++++++++++++ tests/test_simulation_resources.py | 140 +++++++ tests/test_spectral_eqe.py | 359 ++++++++++++++++++ 10 files changed, 1628 insertions(+), 1 deletion(-) create mode 100644 docs/ENSEMBLE_RESOURCE_ESTIMATES.md create mode 100644 docs/SPECTRAL_EQE_AND_DATA.md create mode 100644 simulation_resources.py create mode 100644 spectral_eqe.py create mode 100644 tests/test_simulation_resources.py create mode 100644 tests/test_spectral_eqe.py diff --git a/.github/workflows/plane-wave-preview.yml b/.github/workflows/plane-wave-preview.yml index 8723f7c..7ea44ae 100644 --- a/.github/workflows/plane-wave-preview.yml +++ b/.github/workflows/plane-wave-preview.yml @@ -60,6 +60,7 @@ jobs: cd "$RUNNER_TEMP" "$RUNNER_TEMP/buyer/bin/mongelli-plane-wave" --wavelength-count 3 --angle-count 3 --output preview.png "$RUNNER_TEMP/buyer/bin/python" -c "from dipole_reference import SingleInterface; r=SingleInterface(550, 20, 1.0).decay_rates(); assert r.vertical == r.parallel == 1.0" + "$RUNNER_TEMP/buyer/bin/python" -c "from spectral_eqe import OpticalRateBudget; from simulation_resources import pixel_study; b=OpticalRateBudget([550], [1], [.5], 'homogeneous'); assert b.monochromatic_eqe(intrinsic_yield=1, charge_balance=1, exciton_yield=1)[0] == .5; assert len(pixel_study()['cases']) == 4" test -s preview.png && test -s preview.csv && test -s preview.json - name: Audit installed runtime dependencies and produce CycloneDX inventory run: | diff --git a/MANIFEST.in b/MANIFEST.in index 3c2b2c2..4576631 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1 +1,3 @@ include docs/DIPOLE_FORMULATION.md +include docs/SPECTRAL_EQE_AND_DATA.md +include docs/ENSEMBLE_RESOURCE_ESTIMATES.md diff --git a/README.md b/README.md index 4be20a2..4a36da3 100644 --- a/README.md +++ b/README.md @@ -125,6 +125,21 @@ license presentation, final screenshots, seller verification and publication. ## Functions and engineering applications +### Dispersive data, spectral EQE and ensemble planning + +See [native optical data and spectral EQE derivation](docs/SPECTRAL_EQE_AND_DATA.md) +and [orientation/resource planning](docs/ENSEMBLE_RESOURCE_ESTIMATES.md). +The new APIs provide pinned native refractiveindex.info ingestion, local +measurement CSV ingestion, checksum-checked offline catalog records, +non-extrapolating n+ik evaluation, wavelength-resolved layer stacks, explicit +EQE rate composition, incoherent orientation comparisons, and CPU/GPU/RAM +workload estimates. No cloud account is required or provisioned. + +**Scope boundary:** stack-dependent plane-wave R/T/A is wired; finite-layer +dipole emission/collected-flux physics is not yet implemented. EQE composition +requires externally derived physical source-work/escape-rate budgets, not +plane-wave T. The resource estimator is not an implemented GPU/FDTD solver. + ### Research dipole reference API `dipole_reference` provides a separate homogeneous and single-interface diff --git a/docs/ENSEMBLE_RESOURCE_ESTIMATES.md b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md new file mode 100644 index 0000000..3e09b08 --- /dev/null +++ b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md @@ -0,0 +1,169 @@ +# Orientation ensembles and CPU/GPU/RAM planning + +Prepared 2026-10-16. APIs: `spectral_eqe.py` and `simulation_resources.py`. +Estimates are explicit workload models, NOT a GPU implementation, hardware +performance guarantee, measured peak memory, or full-pixel field simulation. + +## Three distinct source problems + +1. **Incoherent dipoles in one plane:** common emitting depth; unrelated + phases; sum emitted powers, not amplitudes. +2. **Incoherent volume ensemble:** distribution of depths and orientations. + In an infinite planar isotropic stack, lateral translation does not change + the optical response; depth and orientation distributions still matter. +3. **Coherent plane-wave illumination or phased emitter plane:** this is not + the same ensemble. It has different source correlations, energy normalization, + interference and computational requirements. The existing plane-wave preview + computes external illumination, not emission. + +For azimuth-integrated power in planar isotropic media, compute vertical and +parallel basis responses once per wavelength and depth. A tilted dipole has +vertical power weight cos(theta)^2. Random isotropic orientations have expected +vertical weight 1/3 and parallel weight 2/3. Horizontal-only, vertical-only and +tilted distributions can reuse these bases. Direction-resolved collection +in a restricted azimuth, anisotropic media or coherent emitters need more than +this two-basis integrated-power model. + +```python +from spectral_eqe import random_dipoles, dipole_ensemble_eqe + +directions = random_dipoles(10000, seed=42) # uniform solid angle, reproducible +# vertical_budget and horizontal_budget must share wavelength/reference. +comparison = dipole_ensemble_eqe( + vertical_budget, horizontal_budget, directions, + intrinsic_yield=0.8, charge_balance=0.9, exciton_yield=1.0, +) +``` + +The result compares equal-excitation **independent fixed dipoles** (average +their yields) with **rapid reorientation** (average rates then form the yield). +These can differ; an orientation-averaged Purcell factor alone does not imply +the same EQE for both ensembles. Correlations, exciton transfer, saturation, +orientation-dependent intrinsic emission and superradiance are not included. + +For a uniform laterally infinite stack, area increases emitter count, photons +and current together; normalized EQE does not inherently scale with area. +Finite pixel edges, contacts, scattering, current crowding and local heating +break that simplification. No finite-pixel EQE is claimed here. + +## Planar workload model + +Let L be wavelength count, Z sampled depths, C thickness configurations, +J finite layers, Q quadrature evaluations per orientation including refinements. + + basis integrals = 2 * L * Z * C + layer-response evaluations = 2 * L * Z * C * Q * J + +The estimator reports an adjustable arithmetic count, a streamed workspace, +retained rate-array storage, and optional N-dipole direction storage/averaging +cost. The workspace omits Python/SciPy interpreter/allocator overhead and any +large field/angle output. It is not process peak RAM. A streaming analytic +orientation distribution needs no N-by-wavelength matrix. + +Increasing physical area or active-emitter density does not multiply these +planar basis integrals. A Cartesian sweep of thicknesses CAN explode: +C is the product of samples per swept layer, not the sum. Adaptive quadrature +cost depends on optical constants, separations and resonances; Q is not fixed +by geometry alone. Array count/FLOP factors are assumptions, not measured code. + +The implemented single-interface kernels use scalar SciPy adaptive quadrature +on the CPU. No CUDA/ROCm/GPU kernel exists. Independent wavelengths/depths can +eventually be parallelized, but GPU performance cannot be inferred from the +device model name or CPU core count. + +## Full 3-D pixel FDTD planning model + +The alternative estimator sizes a uniform 3-D Yee grid with PML added to every +face. Defaults for comparison: + +- Pixel x/y side: 1, 10, 50, 100 um; simulated z: 1 um. +- Active emitter region: 0.1 um; illustrative active density 1000 per um^3. + This is not a measured exciton density or total molecular density. +- Grid spacing: 10 nm; 10 PML cells per face; single-precision real fields. +- 24 real arrays per cell, plus 30% memory headroom. + Six are electromagnetic field components; auxiliary dispersive/PML storage + is budgeted, not derived for a specific material implementation. +- Physical duration: 200 fs; Courant factor 0.99. +- 20 SAMPLE source positions, 3 independent dipole axes each, run sequentially. + This is NOT every emitter in the volume and NOT a demonstrated converged + sampling density. Finite pixels require converged lateral AND depth sampling. +- 300 assumed FLOPs and 192 bytes of traffic per cell-step. Change these for + the implemented update scheme; double precision also needs a correspondingly + revised traffic assumption, not merely doubled allocated field storage. + +For spacing delta and minimum nondispersive index n_min, the reference +Courant timestep is `0.99*n_min*delta/(c*sqrt(3))`. The fastest medium controls +this bound; a maximum-index CFL formula would be wrong. Dispersive models can +require smaller timesteps. Spatial resolution must also resolve wavelength +inside the largest-index medium, skin depths, thin films and geometric edges. +Ten nm/200 fs is a scenario, not a convergence certification. + +| Pixel side (um) | Active volume (um^3) | Planar sheet basis integrals | Planar 20-depth integrals | FDTD provisioned grid GiB | +|---|---:|---:|---:|---:| +| 1 | 0.1 | 62 | 1,240 | 0.201 | +| 10 | 10 | 62 | 1,240 | 14.511 | +| 50 | 250 | 62 | 1,240 | 351.482 | +| 100 | 1,000 | 62 | 1,240 | 1,400.333 | + +The planar rows assume 31 wavelengths and one thickness configuration. +FDTD RAM is for ONE sequential run, not 60 resident copies. Work multiplies +by independent source runs; memory multiplies only if those runs are simultaneous. +Summing many phase-coherent sources in one run is not an incoherent ensemble. +Stochastic-source/reciprocity methods can reduce work but require separate +implementation and statistical validation. + +The grid estimate excludes full-field snapshots, frequency-domain monitors, +geometry metadata, MPI/halo copies and CPU/GPU transfers. Volume snapshots +can easily dominate storage. A full-spectrum FDTD pulse can cover several +wavelengths in one run if duration/sampling and dispersion are adequate; +do not multiply its cost by planar wavelength count automatically. + +## This workstation and honest runtime limits + +Read-only inspection found an Intel Core i9-13900HX (24 cores, 32 logical +processors), approximately 31.69 GiB system RAM and an RTX 4070 Laptop GPU +reporting 8188 MiB VRAM. These are observed capacity facts, not measured +sustained bandwidth or floating-point throughput. + +Under the table's model, the 10 um pixel exceeds that GPU's memory but fits +the nominal system RAM budget; other application allocations/OS overhead can +still make it impractical. The 50/100 um cases exceed both capacities by large +factors. No multi-GPU capability is implied. + +`Hardware` optionally accepts USER-SUPPLIED sustained CPU/GPU GFLOP/s and GB/s. +The estimator then reports only the idealized lower bound: + + seconds >= max(assumed_FLOPs / sustained_FLOPs_per_second, + assumed_memory_traffic / sustained_bytes_per_second) + +This roofline lower bound is not elapsed-time prediction. If a grid does not +fit, the single-device bound does not mean the run is feasible. Without supplied +rates, no CPU/GPU runtime is invented. The planning tool itself does not execute +FDTD on either processor. + +## Reproduce the scenario report and CPU reference benchmark + +```powershell +python -m simulation_resources --output pixel-resource-study.json --benchmark-cpu +``` + +Or use installed `mongelli-resource-study` with the same flags. Output creation +does not overwrite existing files. `--side-um`, `--emitter-depth-um`, and +`--density-per-um3` customize selected dimensions; use Python APIs for all +quadrature, geometry, precision and hardware assumptions. + +The optional CPU benchmark measures only the implemented single-interface +reference at 550 nm/20 nm for matched and synthetic lossy substrates, with +warmup and five sequential repeats. It reports median/min/max wall time and +quadrature counts. It does NOT benchmark multilayer OLED emission, full +volumetric FDTD, GPU execution, concurrency scaling or peak memory. +It also composes horizontal, vertical, 45-degree and isotropic source-work +responses from these two computed bases, including a 10,000-dipole sample +with seed 42. These are electromagnetic decay/source-work comparisons, +not collected photon fractions or EQE. The lossy reference has relative +source work approximately 2.07018 horizontally, 7.75800 vertically and +3.96612 for an analytic isotropic ensemble. + +For actual procurement/runtime estimates, next benchmark the completed source +solver on representative material stacks and convergence targets, then implement +and measure any GPU backend. Retain workload definitions and uncertainty ranges. diff --git a/docs/SPECTRAL_EQE_AND_DATA.md b/docs/SPECTRAL_EQE_AND_DATA.md new file mode 100644 index 0000000..4a44de8 --- /dev/null +++ b/docs/SPECTRAL_EQE_AND_DATA.md @@ -0,0 +1,210 @@ +# Dispersive stacks, native optical data and spectral EQE + +Research milestone, 2026-10-16. Implemented foundations are in `spectral_eqe.py`. +**There is not yet a finite emitting-layer multilayer source solver.** +The data catalog, stack evaluator, plane-wave checks and EQE rate composition +are working components, not a complete device-performance prediction. + +## What wavelength-specific EQE means + +Let a stack be defined by thicknesses d_j and relative complex indices +N_j(lambda)=n_j(lambda)+i k_j(lambda), under exp(-i omega t). +For a specified emitter position, dipole orientation and collection geometry, +the eventual emitting-source solver must compute: + + F_total(lambda; {d_j,N_j}, z, orientation) = P_EM / P_host + F_escape(lambda; {d_j,N_j}, z, orientation) = P_collected / P_host + +P_EM is source work into all electromagnetic channels, including loss; +P_collected is emitted flux into a specified external collection region. +Both use the SAME lossless homogeneous-host reference for the same dipole. +They are not the R/T coefficients of an externally incident plane wave. +Air escape, substrate emission and detector-NA collection are different regions. + +Let q0 be intrinsic photon yield in that reference environment, eta_charge +the charge-balance factor and eta_exciton the fraction of injected pairs +forming the modeled emissive state. For an independently evaluated +monochromatic transition, with the intrinsic nonradiative rate held fixed: + + EQE_mono(lambda) = + eta_charge * eta_exciton * q0 * F_escape(lambda) + / [1 - q0 + q0 * F_total(lambda)] + +This follows from Gamma_rad,host/(Gamma_rad,host+Gamma_nr)=q0, +Gamma_EM=F_total*Gamma_rad,host and +Gamma_escape=F_escape*Gamma_rad,host. q0=0 yields no photons; q0=1 +reduces to the escape/total optical-rate fraction (unless total rate vanishes, +in which case the steady-state yield is undefined). + +Thickness and n+ik affect both optical rates. They do not determine +charge balance, exciton formation, intrinsic emitter loss or spatial excitation. +No universal 25% spin factor, unit charge balance or unit q0 is silently assumed. + +### A broadband emitter has a different denominator + +For one emitter population with intrinsic PHOTON spectral density S0(lambda) +normalized to integral S0 d(lambda)=1: + + D = 1 - q0 + q0 * integral S0(lambda) F_total(lambda) d(lambda) + d(EQE)/d(lambda) = + eta_charge * eta_exciton * q0 * S0(lambda) F_escape(lambda) / D + +Integrate the latter to get dimensionless broadband EQE. It is incorrect +to apply an independent lifetime denominator to every spectral bin of this +same emitter. A spectrum in optical power per nm must be converted to photon +weights (multiply by lambda/(h*c)) and explicitly normalized first. +The supplied wavelength range must cover the intended emission spectrum; +truncated data normalized to one describe a conditional/truncated model. +The implementation requires a nonnegative already-normalized photon density; +it never silently renormalizes, extrapolates or mixes power/photon units. + +## Native material-data catalog + +`OpticalDataLake(root)` uses an explicit local or already-mounted filesystem +root. It stores immutable dataset manifests plus SHA-256-addressed original +source bytes. Manifest/source checksums are checked on load. These detect +accidental mutation, not malicious alteration of both object and checksum; +they are not cryptographic signatures or proof of material accuracy. + +Implemented providers: + +1. Direct, explicit HTTPS fetch of individual native + [refractiveindex.info database](https://github.com/polyanskiy/refractiveindex.info-database) + YAML datasets pinned to a complete Git commit. +2. Native YAML import from an explicitly identified copy of that database. +3. Local measurements or other licensed sources via exact CSV headers + `wavelength_nm,n,k`, with required source, license, citation and conditions. + +This is a local-first curated catalog, not a provisioned Azure/GCP/S3 service, +full database mirror, or integration with every refractive-index vendor. +New providers can import that CSV schema without losing source provenance. +An existing mounted data lake can be the root; cloud authorization, object-store +adapters, access governance and multi-user transactional writes are future work. +Do not point it at untrusted manifests or overwrite scientific source objects. + +### Supported native formats and units + +- `tabulated nk`, or separate `tabulated n` and `tabulated k`. +- Sellmeier `formula 1` and `formula 2` for n, optionally with tabulated k. +- Unsupported analytic formulas, anisotropic/tensor models, nonlinear n2 + and duplicate components raise errors rather than being approximated. +- Database wavelengths in micrometers are converted to nm. n and k tables + are linearly interpolated separately, only inside their overlapping domains. +- Missing k is NOT evidence of zero absorption. An n-only dataset needs + explicit `assume_lossless=True`, which is preserved in metadata. +- No automatic material selection, zero-loss fallback, extrapolation, or + network access occurs during stack evaluation. + +The upstream database declares CC0-1.0; retain database/original publication +citations and measurement conditions. That does not make its hosted articles, +figures or unrelated commercial databases freely redistributable. Import +only data you have rights to use. Sample preparation, film thickness, +temperature, deposition method, anisotropy and wavelength range affect suitability. + +### Example: pinned measured silver data + +```python +from pathlib import Path +from spectral_eqe import OpticalDataLake + +lake = OpticalDataLake(Path.home() / "optical-material-catalog") +commit = "c5c2f188e848453def5970e347399d653df2ffc2" +silver = lake.fetch_rii( + "Ag-Johnson-1972", + "database/data/main/Ag/nk/Johnson.yml", + commit, +) +print(silver.index([450.0, 550.0, 650.0])) # n + i*k +print(silver.metadata["citation"]) +``` + +The first fetch downloads only this public dataset. Subsequent same-source +calls use verified local files without network access. Reusing an ID for +another revision or lossless assumption is rejected; use a new ID. +Select your own root rather than committing a private material lake. + +Fused silica is also supported through +`database/data/main/SiO2/nk/Malitson.yml`, using formula 1; treating its +missing k as zero requires an explicit lossless assumption. Its room-temperature +optical formula is not a substitute for an OLED organic emitter measurement. + +Database citation: +M. N. Polyanskiy, Scientific Data 11, 94 (2024), +[doi:10.1038/s41597-023-02898-2](https://doi.org/10.1038/s41597-023-02898-2). +Upstream [database notice](https://refractiveindex.info/about), +checked 2026-10-16. + +## Layer stack wiring + +```python +from spectral_eqe import Layer, DispersiveStack + +# Each ID must already be imported from appropriate measurement data. +stack = DispersiveStack( + layers=( + Layer("ambient", None), # semi-infinite top half-space + Layer("transport-top", 35.0), + Layer("emitter-host", 25.0), + Layer("transport-bottom", 40.0), + Layer("metal-electrode", 80.0), + Layer("substrate", None), # semi-infinite bottom half-space + ), + emitter_layer=2, + emitter_depth_nm=12.5, # from the top of that layer +) +evaluated = stack.evaluate(lake, [450.0, 500.0, 550.0, 600.0]) +print(evaluated["indices"]) # [wavelength, layer], complex +print(evaluated["thicknesses_nm"]) + +diagnostic = stack.plane_wave_spectrum(lake, [450.0, 500.0, 550.0, 600.0]) +# diagnostic["R"], ["T"], ["A"] are NOT an EQE result. +``` + +Interior thicknesses must be nonnegative; the emitter position must be strictly +inside a positive-thickness emitting layer. For the current model contract, +the emitting host and exterior half-spaces must be lossless positive-index +media. Absorbing interior non-emitting films are allowed. +Every material's provenance and explicit lossless assumptions accompany +the evaluated stack. Sweeps construct new stacks, not mutate catalog datasets. + +## EQE composition API + +```python +from spectral_eqe import OpticalRateBudget + +# ANALYTIC homogeneous example, not output from a finite OLED device solver. +# Total source work is 1; one homogeneous hemisphere collects 1/2. +budget = OpticalRateBudget( + wavelengths_nm=[500.0, 550.0, 600.0], + total=[1.0, 1.0, 1.0], + escape=[0.5, 0.5, 0.5], + reference="Analytic homogeneous lossless host, one hemisphere, same dipole", +) +print(budget.monochromatic_eqe( + intrinsic_yield=0.8, charge_balance=0.9, exciton_yield=1.0, +)) # 0.36 at each wavelength for this stipulated model +``` + +The constructor checks rate bounds/grid/reference, not whether a caller's +source calculation is scientifically valid. Raw rates MUST be derived and +validated independently; no input label alone proves this. + +`orientation_average` combines total and escape powers before calculating a +rapidly reorienting population's yield. Independently excited, fixed-orientation +dipoles have different lifetimes: compute each one's yield then average +populations. `dipole_ensemble_eqe` returns both results explicitly rather than +silently assuming they are equivalent. It currently treats monochromatic +transition arrays, planar isotropic media and azimuth-integrated collection. +Broadband, spatially distributed fixed-emitter populations require per-population +integrated rate denominators before averaging; they are not implemented here. + +## Remaining finite-stack physics + +The stack is now wired to dispersive data, but it is **not yet wired to a +finite emitting-layer source solver**. Next implement and validate stable +multiple reflection for dipole source terms between both interfaces, +orientation-resolved source work, emitted far-field flux, absorption and +guided-channel accounting with energy/convergence diagnostics. +Only then produce layer-thickness-dependent EQE from a stack automatically. +The existing single-interface source-work kernel alone supplies no collected +power, and plane-wave transmittance must not fill that gap. diff --git a/pyproject.toml b/pyproject.toml index 87e3695..67aa171 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,6 +16,7 @@ dependencies = [ "matplotlib>=3.7,<4", "tmm>=0.2.0,<0.3", "scipy>=1.11,<2", + "PyYAML>=6.0.2,<7", ] classifiers = [ "Development Status :: 3 - Alpha", @@ -28,12 +29,13 @@ dev = ["pytest>=8,<10", "build>=1.2,<2", "twine>=6,<7", "pip-audit>=2.9,<3"] [project.scripts] mongelli-plane-wave = "tmm_OLED9:main" +mongelli-resource-study = "simulation_resources:main" [project.urls] Repository = "https://github.com/gmongell/Mathematica_python_TMM_OLEDPixel" [tool.setuptools] -py-modules = ["tmm_OLED9", "dipole_reference"] +py-modules = ["tmm_OLED9", "dipole_reference", "spectral_eqe", "simulation_resources"] [tool.pytest.ini_options] testpaths = ["tests"] diff --git a/simulation_resources.py b/simulation_resources.py new file mode 100644 index 0000000..99a6009 --- /dev/null +++ b/simulation_resources.py @@ -0,0 +1,266 @@ +"""Transparent workload/RAM estimates, not a GPU solver or measured speed prediction.""" + +from dataclasses import dataclass +import argparse +import json +import math +from numbers import Real +from pathlib import Path +import statistics +import time + + +def _positive(name: str, value: float) -> None: + if not isinstance(value, Real) or not math.isfinite(value) or value <= 0: + raise ValueError(f"{name} must be finite and positive") + + +def _count(name: str, value: int) -> None: + if isinstance(value, bool) or not isinstance(value, int) or value < 1: + raise ValueError(f"{name} must be a positive integer") + + +@dataclass(frozen=True) +class Hardware: + """User-supplied SUSTAINED rates, not vendor peak/nameplate specifications.""" + + cpu_gflops: float + cpu_bandwidth_gb_s: float + gpu_gflops: float + gpu_bandwidth_gb_s: float + system_ram_gib: float + gpu_vram_gib: float + + def __post_init__(self): + for name, value in vars(self).items(): + _positive(name, value) + + +def _roofline(flops: float, traffic_bytes: float, gflops: float, gb_s: float) -> float: + return max(flops / (gflops * 1e9), traffic_bytes / (gb_s * 1e9)) + + +def planar_estimate( + *, wavelengths: int, finite_layers: int, emitter_depths: int, + quadrature_evaluations: int, thickness_configurations: int = 1, + ensemble_emitters: int = 1, bytes_per_real: int = 8, + flops_per_layer_evaluation: int = 200, +) -> dict: + """Two orientation bases per depth; uniform lateral area needs no new solves. + + Quadrature count is per orientation and includes convergence refinements. + Assumes streaming wavelength/depth integration, no full-field spatial output. + """ + arguments = locals().copy() + for name, value in arguments.items(): + _count(name, value) + if bytes_per_real not in (4, 8): + raise ValueError("Real precision must be 4 or 8 bytes") + solves = 2 * wavelengths * emitter_depths * thickness_configurations + layer_evaluations = solves * quadrature_evaluations * finite_layers + flops = layer_evaluations * flops_per_layer_evaluation + # Workspace, index tables and a reusable pair of rate spectra; not per emitter. + working_bytes = bytes_per_real * ( + 16 * (finite_layers + 2) + 2 * wavelengths * (finite_layers + 2) + 8 * wavelengths + ) + retained_bytes = bytes_per_real * 4 * wavelengths * emitter_depths * thickness_configurations + return { + "method": "planar_incoherent_two_orientation_basis", + "basis_integrals": solves, + "layer_evaluations": layer_evaluations, + "assumed_flops": flops, + "streamed_workspace_bytes": working_bytes, + "retained_rate_output_bytes": retained_bytes, + "sampled_orientation_storage_bytes_if_retained": 3 * bytes_per_real * ensemble_emitters, + "orientation_composition_operations": 20 * wavelengths * ensemble_emitters, + "assumptions": arguments, + "runtime_status": "No seconds predicted: complex recursion/quadrature throughput must be benchmarked", + "gpu_status": "No GPU implementation; scalar adaptive SciPy quadrature currently executes on CPU", + } + + +def emitter_population(*, area_um2: float, depth_um: float, density_per_um3: float) -> dict: + for name, value in locals().copy().items(): + _positive(name, value) + volume = area_um2 * depth_um + expected = volume * density_per_um3 + if not math.isfinite(expected): + raise ValueError("Emitter population overflowed") + return {"area_um2": area_um2, "volume_um3": volume, "expected_emitters": expected} + + +def fdtd_estimate( + *, size_um: tuple[float, float, float], cell_nm: float, duration_fs: float, + minimum_index: float = 1.0, pml_cells: int = 10, bytes_per_real: int = 4, + allocated_real_arrays_per_cell: int = 24, flops_per_cell_step: int = 300, + bytes_transferred_per_cell_step: int = 192, independent_source_runs: int = 1, + hardware: Hardware | None = None, +) -> dict: + """Uniform 3-D Yee-grid planning model with configurable ADE/PML array budget. + + CFL uses the minimum nondispersive refractive index (fastest phase velocity), + NOT the largest index. Dispersive media can impose stricter stability bounds. + """ + if len(size_um) != 3: + raise ValueError("Physical domain must have three dimensions in micrometers") + for size in size_um: + _positive("Domain dimension (um)", size) + for name, value in (("Cell size (nm)", cell_nm), ("Duration (fs)", duration_fs), ("Minimum index", minimum_index)): + _positive(name, value) + for name, value in ( + ("PML cells", pml_cells), ("Array count", allocated_real_arrays_per_cell), + ("FLOPs per cell step", flops_per_cell_step), + ("Bytes transferred per cell step", bytes_transferred_per_cell_step), + ("Independent source runs", independent_source_runs), + ): + _count(name, value) + if allocated_real_arrays_per_cell < 6 or bytes_per_real not in (4, 8): + raise ValueError("Need at least six field arrays and 4- or 8-byte real precision") + shape = tuple(math.ceil(size * 1000 / cell_nm) + 2 * pml_cells for size in size_um) + cells = math.prod(shape) + dt_fs = 0.99 * minimum_index * cell_nm * 1e-9 / (299792458 * math.sqrt(3)) * 1e15 + steps = math.ceil(duration_fs / dt_fs) + updates = cells * steps * independent_source_runs + raw_bytes = cells * allocated_real_arrays_per_cell * bytes_per_real + provision_bytes = math.ceil(raw_bytes * 1.3) + result = { + "method": "uniform_3d_fdtd_planning_only", + "grid_shape_including_pml": shape, "cells": cells, + "courant_step_fs": dt_fs, "time_steps_per_run": steps, + "cell_updates_all_runs": updates, + "raw_grid_bytes": raw_bytes, "provisioned_grid_bytes_30pct_headroom": provision_bytes, + "assumed_flops": updates * flops_per_cell_step, + "assumed_memory_traffic_bytes": updates * bytes_transferred_per_cell_step, + "independent_source_runs": independent_source_runs, + "limitations": ( + "No FDTD implementation. Excludes full-field snapshots, DFT monitors, geometry, " + "MPI copies, device transfers, dispersive stability analysis and convergence runs. " + "Array/FLOP/traffic counts are user-adjustable planning assumptions. " + "One coherent source run does not represent an incoherent emitter ensemble." + ), + } + if hardware is not None: + result.update({ + "cpu_ideal_roofline_seconds_lower_bound": _roofline(result["assumed_flops"], result["assumed_memory_traffic_bytes"], hardware.cpu_gflops, hardware.cpu_bandwidth_gb_s), + "gpu_ideal_roofline_seconds_lower_bound": _roofline(result["assumed_flops"], result["assumed_memory_traffic_bytes"], hardware.gpu_gflops, hardware.gpu_bandwidth_gb_s), + "fits_cpu_ram_budget": provision_bytes <= hardware.system_ram_gib * 2**30, + "fits_gpu_vram_budget": provision_bytes <= hardware.gpu_vram_gib * 2**30, + "hardware_assumptions": vars(hardware), + }) + return result + + +def pixel_study( + *, side_lengths_um: tuple[float, ...] = (1, 10, 50, 100), + emitter_depth_um: float = .1, active_density_per_um3: float = 1000, + simulation_depth_um: float = 1, volume_depth_samples: int = 20, + fdtd_source_positions: int = 20, hardware: Hardware | None = None, +) -> dict: + """Compare sheet/volume planar ensembles with a finite-pixel reference workload.""" + _count("Volume depth samples", volume_depth_samples) + _count("FDTD sampled source positions", fdtd_source_positions) + _positive("Simulation depth (um)", simulation_depth_um) + if not side_lengths_um: + raise ValueError("At least one pixel side length is required") + if emitter_depth_um > simulation_depth_um: + raise ValueError("Emission region cannot exceed the simulated domain depth") + rows = [] + for side in side_lengths_um: + _positive("Pixel side (um)", side) + population = emitter_population( + area_um2=side**2, depth_um=emitter_depth_um, density_per_um3=active_density_per_um3, + ) + common = dict(wavelengths=31, finite_layers=6, quadrature_evaluations=2000, + ensemble_emitters=max(1, math.ceil(population["expected_emitters"]))) + rows.append({ + "pixel_side_um": side, **population, + "emitter_sheet": planar_estimate(**common, emitter_depths=1), + "emitter_volume": planar_estimate(**common, emitter_depths=volume_depth_samples), + "finite_pixel_fdtd": fdtd_estimate( + size_um=(side, side, simulation_depth_um), cell_nm=10, duration_fs=200, + independent_source_runs=3 * fdtd_source_positions, hardware=hardware, + ), + }) + return { + "status": "Planning estimates, not an implemented finite-stack or full-wave simulation", + "assumptions": { + "planar_wavelengths": 31, "finite_layers": 6, + "quadrature_evaluations_per_orientation": 2000, + "volume_depth_samples": volume_depth_samples, + "active_density_per_um3": active_density_per_um3, + "fdtd_source_positions": fdtd_source_positions, + "fdtd_source_runs_per_position": 3, + "fdtd_cell_nm": 10, "fdtd_duration_fs": 200, + "fdtd_simulation_depth_um": simulation_depth_um, + "fdtd_real_arrays_per_cell": 24, "fdtd_precision_bytes": 4, + "fdtd_pml_cells_per_face": 10, + }, + "cases": rows, + } + + +def benchmark_cpu_reference(*, repeats: int = 5) -> dict: + """Measure only the implemented single-interface CPU kernels on this machine.""" + from dipole_reference import SingleInterface + from spectral_eqe import random_dipoles + + _count("Benchmark repeats", repeats) + directions = random_dipoles(10000, seed=42) + sampled_vertical_fraction = float((directions[:, 2] ** 2).mean()) + records = [] + for name, epsilon in [("matched_homogeneous", 1.0), ("synthetic_lossy_interface", -5 + .5j)]: + interface = SingleInterface(550, 20, epsilon) + interface.decay_rates() + samples = [] + for _ in range(repeats): + start = time.perf_counter() + result = interface.decay_rates() + samples.append(time.perf_counter() - start) + records.append({ + "case": name, "wavelength_nm": 550, "distance_nm": 20, + "median_wall_seconds": statistics.median(samples), + "min_wall_seconds": min(samples), "max_wall_seconds": max(samples), + "repeats": repeats, "vertical_source_work": result.vertical, + "parallel_source_work": result.parallel, + "quadrature_evaluations_per_call": result.evaluations, + "orientation_source_work_not_eqe": { + "vertical": result.vertical, + "horizontal_plane": result.parallel, + "tilted_45_degrees": .5 * result.vertical + .5 * result.parallel, + "isotropic_analytic": result.vertical / 3 + 2 * result.parallel / 3, + "isotropic_sampled": sampled_vertical_fraction * result.vertical + (1 - sampled_vertical_fraction) * result.parallel, + }, + }) + return { + "method": "Measured sequential CPU single-interface reference; no finite-stack extrapolation", + "cases": records, + "orientation_sampling": { + "count": len(directions), "seed": 42, + "mean_vertical_fraction": sampled_vertical_fraction, + "quantity": "Total electromagnetic source work / homogeneous-host power, not collected flux or EQE", + }, + "gpu": "Not measured: no GPU kernel is implemented", + "memory": "No process peak-RAM measurement; use the separately labeled planning model", + } + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--side-um", nargs="+", type=float, default=[1, 10, 50, 100]) + parser.add_argument("--emitter-depth-um", type=float, default=.1) + parser.add_argument("--density-per-um3", type=float, default=1000) + parser.add_argument("--benchmark-cpu", action="store_true") + args = parser.parse_args() + report = pixel_study(side_lengths_um=tuple(args.side_um), emitter_depth_um=args.emitter_depth_um, + active_density_per_um3=args.density_per_um3) + if args.benchmark_cpu: + report["measured_reference_cpu"] = benchmark_cpu_reference() + with args.output.open("x", encoding="utf-8") as output: + json.dump(report, output, indent=2, allow_nan=False) + output.write("\n") + print(f"Created planning report: {args.output}; no GPU execution or runtime measurement claimed") + + +if __name__ == "__main__": + main() diff --git a/spectral_eqe.py b/spectral_eqe.py new file mode 100644 index 0000000..6cae336 --- /dev/null +++ b/spectral_eqe.py @@ -0,0 +1,463 @@ +"""Dispersive optical-data catalog and rate-model EQE; no finite-stack emitter solver.""" + +import csv +from dataclasses import dataclass +import hashlib +import io +import json +from pathlib import Path, PurePosixPath +import re +from typing import Mapping +from urllib.request import urlopen + +import numpy as np +from numpy.typing import ArrayLike, NDArray +from scipy.integrate import trapezoid +from tmm import coh_tmm +import yaml + + +def _grid(values: ArrayLike) -> NDArray[np.float64]: + grid = np.asarray(values, dtype=float) + if grid.ndim != 1 or grid.size == 0 or not np.all(np.isfinite(grid)) or np.any(grid <= 0): + raise ValueError("Wavelengths must be a nonempty finite positive 1-D array (nm)") + if np.any(np.diff(grid) <= 0): + raise ValueError("Wavelengths must be strictly increasing") + return grid + + +def _probability(name: str, value: float) -> None: + if not np.isfinite(value) or not 0 <= value <= 1: + raise ValueError(f"{name} must be finite and in [0, 1]") + + +def _table(values: ArrayLike, columns: int) -> NDArray[np.float64]: + data = np.asarray(values, dtype=float) + if data.ndim != 2 or data.shape[1] != columns or len(data) < 2 or not np.all(np.isfinite(data)): + raise ValueError("Optical table must have at least two finite rows and the required columns") + _grid(data[:, 0]) + if np.any(data[:, 1:] < 0): + raise ValueError("Passive optical tables require n >= 0 and k >= 0") + return data + + +def _within(grid: NDArray[np.float64], lower: float, upper: float) -> None: + if np.any(grid < lower) or np.any(grid > upper): + raise ValueError(f"Wavelength outside dataset domain [{lower:g}, {upper:g}] nm; extrapolation forbidden") + + +@dataclass(frozen=True) +class Material: + """A verified catalog record. Data are interpreted only through index().""" + + dataset_id: str + metadata: Mapping + blocks: tuple[dict, ...] + assume_lossless: bool + + def index(self, wavelengths_nm: ArrayLike) -> NDArray[np.complex128]: + grid = _grid(wavelengths_nm) + n = k = None + for block in self.blocks: + kind = block["type"] + if kind in ("tabulated nk", "tabulated n", "tabulated k"): + table = _table(block["rows_nm"], 3 if kind == "tabulated nk" else 2) + _within(grid, table[0, 0], table[-1, 0]) + if kind != "tabulated k": + if n is not None: + raise ValueError("Multiple refractive-index components in dataset") + n = np.interp(grid, table[:, 0], table[:, 1]) + if kind != "tabulated n": + if k is not None: + raise ValueError("Multiple extinction components in dataset") + k = np.interp(grid, table[:, 0], table[:, -1]) + elif kind in ("formula 1", "formula 2"): + if n is not None: + raise ValueError("Multiple refractive-index components in dataset") + lower, upper = np.asarray(block["range_nm"], dtype=float) + if not np.isfinite([lower, upper]).all() or not 0 < lower < upper: + raise ValueError("Formula range must be finite, positive and increasing") + _within(grid, lower, upper) + coefficients = np.asarray(block["coefficients"], dtype=float) + if coefficients.ndim != 1 or len(coefficients) < 3 or len(coefficients) % 2 != 1 or not np.isfinite(coefficients).all(): + raise ValueError("Sellmeier coefficients must be a finite constant plus coefficient pairs") + wavelength_squared = (grid / 1000) ** 2 + n_squared = np.full_like(grid, 1 + coefficients[0]) + for strength, pole in zip(coefficients[1::2], coefficients[2::2]): + denominator = wavelength_squared - (pole**2 if kind == "formula 1" else pole) + if np.any(denominator == 0): + raise ValueError("Wavelength coincides with a dispersion-formula pole") + n_squared += strength * wavelength_squared / denominator + if np.any(n_squared <= 0) or not np.isfinite(n_squared).all(): + raise ValueError("Dispersion formula produced invalid dielectric refractive index") + n = np.sqrt(n_squared) + else: + raise ValueError(f"Unsupported optical-data block: {kind}") + if n is None: + raise ValueError("Dataset has no refractive-index component") + if k is None: + if not self.assume_lossless: + raise ValueError("Extinction coefficient is missing; explicit assume_lossless=True is required") + k = np.zeros_like(n) + result = n + 1j * k + if not np.isfinite(result).all(): + raise ValueError("Material evaluation produced nonfinite optical constants") + return result + + +class OpticalDataLake: + """Immutable content-addressed local catalog; simulations never fetch implicitly. + + The root may be a local directory or an already-mounted shared filesystem. + This is not a cloud provisioning client or a database-wide mirror. + """ + + def __init__(self, root: str | Path): + self.root = Path(root) + self.root.mkdir(parents=True, exist_ok=True) + + @staticmethod + def _identifier(dataset_id: str) -> str: + if not re.fullmatch(r"[A-Za-z0-9][A-Za-z0-9_.-]{0,127}", dataset_id): + raise ValueError("Dataset ID must be a simple catalog identifier, not a path") + return dataset_id + + def list_datasets(self) -> list[str]: + return sorted(path.stem for path in self.root.glob("*.json")) + + def _store(self, dataset_id: str, raw: bytes, blocks: list[dict], metadata: dict, assume_lossless: bool) -> Material: + dataset_id = self._identifier(dataset_id) + for field in ("source", "license", "citation"): + if not isinstance(metadata.get(field), str) or not metadata[field].strip(): + raise ValueError(f"Material provenance requires {field}") + if not isinstance(assume_lossless, bool): + raise ValueError("assume_lossless must be explicitly boolean") + if not blocks: + raise ValueError("Dataset contains no supported optical blocks") + for block in blocks: + if block["type"].startswith("tabulated "): + _table(block["rows_nm"], 3 if block["type"] == "tabulated nk" else 2) + elif block["type"] not in ("formula 1", "formula 2"): + raise ValueError(f"Unsupported optical-data block: {block['type']}") + raw_sha = hashlib.sha256(raw).hexdigest() + record = { + "schema_version": 1, "dataset_id": dataset_id, + "metadata": metadata, "blocks": blocks, + "assume_lossless": assume_lossless, "raw_sha256": raw_sha, + } + payload = json.dumps(record, sort_keys=True, allow_nan=False).encode("utf-8") + envelope = {"payload": record, "sha256": hashlib.sha256(payload).hexdigest()} + path = self.root / f"{dataset_id}.json" + encoded = json.dumps(envelope, indent=2, sort_keys=True, allow_nan=False).encode("utf-8") + if path.exists(): + if path.read_bytes() != encoded: + raise FileExistsError(f"Dataset ID already exists with different content: {dataset_id}") + return self.load(dataset_id) + raw_path = self.root / f"{raw_sha}.source" + if raw_path.exists(): + if hashlib.sha256(raw_path.read_bytes()).hexdigest() != raw_sha: + raise ValueError("Existing raw material object failed checksum validation") + else: + with raw_path.open("xb") as output: + output.write(raw) + with path.open("xb") as output: + output.write(encoded) + return self.load(dataset_id) + + def load(self, dataset_id: str) -> Material: + dataset_id = self._identifier(dataset_id) + envelope = json.loads((self.root / f"{dataset_id}.json").read_text(encoding="utf-8")) + record = envelope["payload"] + payload = json.dumps(record, sort_keys=True, allow_nan=False).encode("utf-8") + if hashlib.sha256(payload).hexdigest() != envelope["sha256"]: + raise ValueError(f"Catalog checksum mismatch: {dataset_id}") + if record["schema_version"] != 1 or record["dataset_id"] != dataset_id: + raise ValueError("Unsupported or mismatched material record") + digest = record["raw_sha256"] + if not re.fullmatch(r"[0-9a-f]{64}", digest): + raise ValueError("Invalid raw material digest") + if hashlib.sha256((self.root / f"{digest}.source").read_bytes()).hexdigest() != digest: + raise ValueError(f"Raw material checksum mismatch: {dataset_id}") + return Material(dataset_id, record["metadata"], tuple(record["blocks"]), record["assume_lossless"]) + + def import_csv( + self, dataset_id: str, path: str | Path, *, source: str, license: str, + citation: str, conditions: str, + ) -> Material: + """Import exact wavelength_nm,n,k headers, preserving original bytes.""" + raw = Path(path).read_bytes() + reader = csv.DictReader(io.StringIO(raw.decode("utf-8-sig"))) + if reader.fieldnames != ["wavelength_nm", "n", "k"]: + raise ValueError("CSV columns must be exactly wavelength_nm,n,k") + rows = [[float(row[key]) for key in reader.fieldnames] for row in reader] + return self._store(dataset_id, raw, [{"type": "tabulated nk", "rows_nm": rows}], { + "source": source, "license": license, "citation": citation, "conditions": conditions, + "format": "CSV wavelength_nm,n,k", + }, False) + + def import_rii( + self, dataset_id: str, raw: bytes, *, source: str, revision: str, + assume_lossless: bool = False, + ) -> Material: + """Import native refractiveindex.info YAML from an explicit source/revision.""" + if not source.strip() or not revision.strip(): + raise ValueError("Source and revision are required") + document = yaml.safe_load(raw) + if not isinstance(document, dict) or not isinstance(document.get("DATA"), list): + raise ValueError("Expected native refractiveindex.info DATA blocks") + blocks = [] + for item in document["DATA"]: + kind = item["type"] + if kind in ("tabulated nk", "tabulated n", "tabulated k"): + columns = 3 if kind == "tabulated nk" else 2 + rows = np.loadtxt(io.StringIO(item["data"]), ndmin=2) + rows = _table(rows, columns).copy() + rows[:, 0] *= 1000 # Database wavelengths are micrometers, API wavelengths are nm. + blocks.append({"type": kind, "rows_nm": rows.tolist()}) + elif kind in ("formula 1", "formula 2"): + blocks.append({ + "type": kind, + "range_nm": [1000 * float(value) for value in str(item["wavelength_range"]).split()], + "coefficients": [float(value) for value in str(item["coefficients"]).split()], + }) + else: + raise ValueError(f"Unsupported refractiveindex.info block: {kind}") + return self._store(dataset_id, raw, blocks, { + "source": source, "revision": revision, "license": "CC0-1.0", + "citation": document.get("REFERENCES", ""), + "conditions": document.get("CONDITIONS", {}), + "comments": document.get("COMMENTS", ""), + "database_citation": "M. N. Polyanskiy, Scientific Data 11, 94 (2024), doi:10.1038/s41597-023-02898-2", + }, assume_lossless) + + def fetch_rii( + self, dataset_id: str, database_path: str, commit: str, *, + assume_lossless: bool = False, + ) -> Material: + """Explicit HTTPS fetch from the official database, pinned to a full Git commit.""" + self._identifier(dataset_id) + if not re.fullmatch(r"[0-9a-f]{40}", commit): + raise ValueError("Use an immutable full lowercase 40-character Git commit") + path = PurePosixPath(database_path) + if (not database_path.startswith("database/data/") or path.suffix != ".yml" + or ".." in path.parts or "\\" in database_path + or not re.fullmatch(r"[A-Za-z0-9_./-]+", database_path)): + raise ValueError("Expected a relative database/data/.../*.yml source path") + url = f"https://raw.githubusercontent.com/polyanskiy/refractiveindex.info-database/{commit}/{database_path}" + if (self.root / f"{dataset_id}.json").exists(): + material = self.load(dataset_id) + if material.metadata["source"] != url or material.assume_lossless != assume_lossless: + raise FileExistsError("Dataset ID is bound to another source or lossless assumption") + return material + with urlopen(url, timeout=30) as response: + raw = response.read(4_000_001) + if len(raw) > 4_000_000: + raise ValueError("Material YAML exceeds the 4 MB per-dataset limit") + return self.import_rii(dataset_id, raw, source=url, revision=commit, assume_lossless=assume_lossless) + + +@dataclass(frozen=True) +class Layer: + material_id: str + thickness_nm: float | None + + +@dataclass(frozen=True) +class DispersiveStack: + """Ordered top-to-bottom stack with an explicit finite lossless emitter layer.""" + + layers: tuple[Layer, ...] + emitter_layer: int + emitter_depth_nm: float + + def __post_init__(self): + if len(self.layers) < 3 or self.layers[0].thickness_nm is not None or self.layers[-1].thickness_nm is not None: + raise ValueError("Stack needs two explicit half-spaces and at least one finite layer") + for layer in self.layers[1:-1]: + if layer.thickness_nm is None or not np.isfinite(layer.thickness_nm) or layer.thickness_nm < 0: + raise ValueError("Interior layer thicknesses must be finite and nonnegative (nm)") + if isinstance(self.emitter_layer, bool) or not isinstance(self.emitter_layer, int) or not 0 < self.emitter_layer < len(self.layers) - 1: + raise ValueError("Emitter must occupy an interior layer") + thickness = self.layers[self.emitter_layer].thickness_nm + if not np.isfinite(self.emitter_depth_nm) or not 0 < self.emitter_depth_nm < thickness: + raise ValueError("Emitter depth must be strictly inside its layer (nm)") + + def evaluate(self, lake: OpticalDataLake, wavelengths_nm: ArrayLike) -> dict: + grid = _grid(wavelengths_nm) + materials = [lake.load(layer.material_id) for layer in self.layers] + indices = np.column_stack([material.index(grid) for material in materials]) + if np.any(indices[:, self.emitter_layer].imag != 0) or np.any(indices[:, self.emitter_layer].real <= 0): + raise ValueError("Emitting host must be lossless and positive index") + if np.any(indices[:, [0, -1]].imag != 0) or np.any(indices[:, [0, -1]].real <= 0): + raise ValueError("Exterior half-spaces must be lossless and positive index") + return { + "wavelengths_nm": grid, + "indices": indices, + "thicknesses_nm": np.array([np.inf if layer.thickness_nm is None else layer.thickness_nm for layer in self.layers]), + "emitter_layer": self.emitter_layer, + "emitter_depth_nm": self.emitter_depth_nm, + "materials": [{"dataset_id": material.dataset_id, **dict(material.metadata), + "assume_lossless": material.assume_lossless} for material in materials], + } + + def plane_wave_spectrum( + self, lake: OpticalDataLake, wavelengths_nm: ArrayLike, *, angle_deg: float = 0, + ) -> dict: + """Bounded dispersive R/T/A integration check; NEVER an emitter EQE result.""" + if not np.isfinite(angle_deg) or not 0 <= angle_deg < 90: + raise ValueError("Incidence angle must be finite and in [0, 90) degrees") + evaluated = self.evaluate(lake, wavelengths_nm) + reflection, transmission = [], [] + for wavelength, indices in zip(evaluated["wavelengths_nm"], evaluated["indices"]): + results = [coh_tmm(p, indices, evaluated["thicknesses_nm"], np.deg2rad(angle_deg), wavelength) + for p in ("s", "p")] + reflection.append(sum(result["R"] for result in results) / 2) + transmission.append(sum(result["T"] for result in results) / 2) + r, t = np.array(reflection), np.array(transmission) + a = 1 - r - t + if not np.isfinite([r, t, a]).all() or np.any(np.array([r, t, a]) < -1e-9): + raise ArithmeticError("Dispersive plane-wave power balance is invalid") + return {**evaluated, "quantity": "plane_wave_power_not_EQE", "R": r, "T": t, "A": a} + + +@dataclass(frozen=True) +class OpticalRateBudget: + """Source-work F_total and collected emitted power F_escape, both / P_host. + + Arrays must come from a physically validated emitting-source calculation, + not from incident-plane-wave transmission. No such finite-stack solver + is manufactured by this data container. + """ + + wavelengths_nm: ArrayLike + total: ArrayLike + escape: ArrayLike + reference: str + + def __post_init__(self): + grid = _grid(self.wavelengths_nm).copy() + total, escape = np.asarray(self.total, dtype=float), np.asarray(self.escape, dtype=float) + if total.shape != grid.shape or escape.shape != grid.shape or not np.isfinite([total, escape]).all(): + raise ValueError("Optical rates must be finite arrays on the wavelength grid") + if np.any(total < 0) or np.any(escape < 0) or np.any(escape > total): + raise ValueError("Require 0 <= escape source-work rate <= total rate") + if not isinstance(self.reference, str) or not self.reference.strip(): + raise ValueError("Optical rate normalization/source reference is required") + total, escape = total.copy(), escape.copy() + for array in (grid, total, escape): + array.setflags(write=False) + object.__setattr__(self, "wavelengths_nm", grid) + object.__setattr__(self, "total", total) + object.__setattr__(self, "escape", escape) + + @classmethod + def orientation_average( + cls, vertical: "OpticalRateBudget", horizontal: "OpticalRateBudget", vertical_fraction: float, + ) -> "OpticalRateBudget": + # This averages rates, not the quantum yields of independently decaying fixed dipoles. + _probability("Vertical fraction", vertical_fraction) + if not np.array_equal(vertical.wavelengths_nm, horizontal.wavelengths_nm): + raise ValueError("Orientation budgets must use exactly the same wavelength grid") + if vertical.reference != horizontal.reference: + raise ValueError("Orientation budgets require the same homogeneous reference normalization") + weight = vertical_fraction + return cls( + vertical.wavelengths_nm, + weight * vertical.total + (1 - weight) * horizontal.total, + weight * vertical.escape + (1 - weight) * horizontal.escape, + vertical.reference, + ) + + def monochromatic_eqe( + self, *, intrinsic_yield: float, charge_balance: float, exciton_yield: float, + ) -> NDArray[np.float64]: + """Independent monochromatic transitions; NOT a broadband spectral density.""" + for name, value in (("Intrinsic yield", intrinsic_yield), ("Charge balance", charge_balance), ("Exciton yield", exciton_yield)): + _probability(name, value) + if intrinsic_yield == 0: + return np.zeros_like(self.total) + denominator = (1 - intrinsic_yield) + intrinsic_yield * self.total + if np.any(denominator <= 0): + raise ValueError("Zero total decay rate: steady-state photon yield is undefined") + return charge_balance * exciton_yield * intrinsic_yield * self.escape / denominator + + def broadband_eqe( + self, intrinsic_photon_density_per_nm: ArrayLike, *, intrinsic_yield: float, + charge_balance: float, exciton_yield: float, + ) -> dict: + """One emitter population: integrate rates BEFORE constructing spectral EQE.""" + for name, value in (("Intrinsic yield", intrinsic_yield), ("Charge balance", charge_balance), ("Exciton yield", exciton_yield)): + _probability(name, value) + spectrum = np.asarray(intrinsic_photon_density_per_nm, dtype=float) + grid = self.wavelengths_nm + if len(grid) < 2 or spectrum.shape != grid.shape or not np.isfinite(spectrum).all() or np.any(spectrum < 0): + raise ValueError("Photon spectrum must be nonnegative and finite on at least two wavelengths") + with np.errstate(over="raise", invalid="raise", divide="raise"): + try: + normalization = trapezoid(spectrum, grid) + except FloatingPointError as error: + raise ValueError("Photon-spectrum integral overflowed") from error + if not np.isclose(normalization, 1, rtol=1e-6, atol=1e-9): + raise ValueError("Intrinsic PHOTON density must integrate to 1 per nm; normalization is explicit") + with np.errstate(over="raise", invalid="raise", divide="raise"): + try: + if intrinsic_yield == 0: + denominator, density = 1.0, np.zeros_like(spectrum) + else: + denominator = (1 - intrinsic_yield) + intrinsic_yield * trapezoid(spectrum * self.total, grid) + if denominator <= 0: + raise ValueError("Zero integrated decay rate: spectral EQE is undefined") + density = charge_balance * exciton_yield * intrinsic_yield * spectrum * self.escape / denominator + integrated_eqe = float(trapezoid(density, grid)) + except FloatingPointError as error: + raise ValueError("Spectral optical-rate integration overflowed") from error + if not np.isfinite(denominator) or not np.isfinite(density).all() or not np.isfinite(integrated_eqe): + raise ValueError("Nonfinite spectral EQE calculation") + return { + "wavelengths_nm": grid.copy(), "eqe_density_per_nm": density, + "integrated_eqe": integrated_eqe, + "rate_denominator": float(denominator), "reference": self.reference, + } + + +def random_dipoles(count: int, *, seed: int) -> NDArray[np.float64]: + """Reproducible uniform orientations on the unit sphere, not coherent phases.""" + if isinstance(count, bool) or not isinstance(count, int) or count < 1: + raise ValueError("Dipole count must be a positive integer") + rng = np.random.default_rng(seed) + cosine = rng.uniform(-1, 1, count) + azimuth = rng.uniform(0, 2 * np.pi, count) + sine = np.sqrt(1 - cosine**2) + return np.column_stack((sine * np.cos(azimuth), sine * np.sin(azimuth), cosine)) + + +def dipole_ensemble_eqe( + vertical: OpticalRateBudget, horizontal: OpticalRateBudget, directions: ArrayLike, *, + intrinsic_yield: float, charge_balance: float, exciton_yield: float, +) -> dict: + """Compare independent fixed dipoles with a rapidly reorienting population. + + Planar isotropic media and azimuth-integrated collection only. Incoherent, + equal excitation probability per emitter; no superradiance or phase array. + Streams emitters to avoid allocating emitter-by-wavelength arrays. + """ + vectors = np.asarray(directions, dtype=float) + if vectors.ndim != 2 or vectors.shape[1] != 3 or len(vectors) == 0 or not np.isfinite(vectors).all(): + raise ValueError("Directions must be a nonempty finite N-by-3 array") + norms = np.linalg.norm(vectors, axis=1) + if np.any(norms == 0) or not np.isfinite(norms).all(): + raise ValueError("Dipole orientation vectors must have finite nonzero norm") + weights = (vectors[:, 2] / norms) ** 2 + averaged = OpticalRateBudget.orientation_average(vertical, horizontal, float(np.mean(weights))) + kwargs = dict(intrinsic_yield=intrinsic_yield, charge_balance=charge_balance, exciton_yield=exciton_yield) + independent = np.zeros_like(averaged.total) + for weight in weights: + budget = OpticalRateBudget.orientation_average(vertical, horizontal, float(weight)) + independent += budget.monochromatic_eqe(**kwargs) / len(weights) + return { + "independent_fixed_dipoles_eqe": independent, + "rapid_reorientation_eqe": averaged.monochromatic_eqe(**kwargs), + "mean_vertical_fraction": float(np.mean(weights)), + "dipole_count": len(weights), + "assumptions": "Incoherent equal-excitation ensemble; planar isotropic media; azimuth-integrated collection", + } diff --git a/tests/test_simulation_resources.py b/tests/test_simulation_resources.py new file mode 100644 index 0000000..f89f71a --- /dev/null +++ b/tests/test_simulation_resources.py @@ -0,0 +1,140 @@ +import math +import json +import sys + +import pytest + +from simulation_resources import Hardware, emitter_population, fdtd_estimate, planar_estimate, pixel_study, benchmark_cpu_reference, main + + +def test_planar_two_bases_reuse_random_ensemble(): + kwargs = dict(wavelengths=31, finite_layers=6, emitter_depths=10, quadrature_evaluations=1000) + one = planar_estimate(**kwargs) + many = planar_estimate(**kwargs, ensemble_emitters=1000000) + assert one["basis_integrals"] == 620 + assert one["layer_evaluations"] == 620 * 1000 * 6 + assert many["assumed_flops"] == one["assumed_flops"] + assert many["streamed_workspace_bytes"] == one["streamed_workspace_bytes"] + assert many["sampled_orientation_storage_bytes_if_retained"] == 24_000_000 + assert many["orientation_composition_operations"] == 20 * 31 * 1000000 + sweep = planar_estimate(**kwargs, thickness_configurations=5) + assert sweep["assumed_flops"] == 5 * one["assumed_flops"] + assert sweep["retained_rate_output_bytes"] == 5 * one["retained_rate_output_bytes"] + + +def test_area_volume_population_scaling(): + small = emitter_population(area_um2=100, depth_um=.1, density_per_um3=1000) + large = emitter_population(area_um2=400, depth_um=.2, density_per_um3=1000) + assert small == {"area_um2": 100, "volume_um3": 10, "expected_emitters": 10000} + assert large["expected_emitters"] == 8 * small["expected_emitters"] + + +def test_pixel_study_compares_area_volume_and_independent_source_work(): + result = pixel_study(side_lengths_um=(1, 10), volume_depth_samples=5, fdtd_source_positions=7) + small, large = result["cases"] + assert large["area_um2"] == 100 * small["area_um2"] + assert large["volume_um3"] == 100 * small["volume_um3"] + assert large["expected_emitters"] == 100 * small["expected_emitters"] + assert small["emitter_volume"]["basis_integrals"] == 5 * small["emitter_sheet"]["basis_integrals"] + assert large["emitter_volume"]["basis_integrals"] == small["emitter_volume"]["basis_integrals"] + assert small["finite_pixel_fdtd"]["independent_source_runs"] == 21 + assert large["finite_pixel_fdtd"]["cells"] > small["finite_pixel_fdtd"]["cells"] + assert "cpu_ideal_roofline_seconds_lower_bound" not in small["finite_pixel_fdtd"] + + +def test_cpu_benchmark_measures_reference_not_full_solver(): + report = benchmark_cpu_reference(repeats=1) + homogeneous, lossy = report["cases"] + assert homogeneous["vertical_source_work"] == homogeneous["parallel_source_work"] == 1 + assert lossy["vertical_source_work"] > 1 + assert lossy["parallel_source_work"] > 0 + assert homogeneous["repeats"] == 1 + assert homogeneous["median_wall_seconds"] > 0 + assert set(homogeneous["orientation_source_work_not_eqe"].values()) == {1.0} + orientations = lossy["orientation_source_work_not_eqe"] + assert orientations["isotropic_analytic"] == pytest.approx(3.966116501134272) + assert orientations["tilted_45_degrees"] == pytest.approx(4.9140864628373) + assert orientations["isotropic_sampled"] == pytest.approx(orientations["isotropic_analytic"], abs=.03) + assert report["orientation_sampling"]["count"] == 10000 + assert report["orientation_sampling"]["seed"] == 42 + assert "no GPU kernel" in report["gpu"] + with pytest.raises(ValueError): + benchmark_cpu_reference(repeats=0) + + +def test_resource_cli_exports_requested_sizes_without_overwrite(tmp_path, monkeypatch): + output = tmp_path / "report.json" + monkeypatch.setattr(sys, "argv", ["resource-study", "--output", str(output), "--side-um", "2", "4"]) + main() + report = json.loads(output.read_text()) + assert [case["pixel_side_um"] for case in report["cases"]] == [2, 4] + assert report["cases"][0]["volume_um3"] == pytest.approx(.4) + previous = output.read_bytes() + with pytest.raises(FileExistsError): + main() + assert output.read_bytes() == previous + + +def test_fdtd_grid_cfl_memory_and_roofline_arithmetic(): + hardware = Hardware(100, 20, 1000, 100, 16, 8) + result = fdtd_estimate(size_um=(1, 2, 3), cell_nm=100, duration_fs=10, + pml_cells=2, hardware=hardware) + assert result["grid_shape_including_pml"] == (14, 24, 34) + cells = 14 * 24 * 34 + assert result["cells"] == cells + assert result["raw_grid_bytes"] == cells * 24 * 4 + assert result["provisioned_grid_bytes_30pct_headroom"] == math.ceil(cells * 24 * 4 * 1.3) + dt = .99 * 100e-9 / (299792458 * math.sqrt(3)) * 1e15 + assert result["courant_step_fs"] == pytest.approx(dt) + assert result["time_steps_per_run"] == math.ceil(10 / dt) + work = result["cell_updates_all_runs"] + assert result["cpu_ideal_roofline_seconds_lower_bound"] == pytest.approx(max(work * 300 / 100e9, work * 192 / 20e9)) + assert result["gpu_ideal_roofline_seconds_lower_bound"] == pytest.approx(max(work * 300 / 1000e9, work * 192 / 100e9)) + assert result["fits_cpu_ram_budget"] and result["fits_gpu_vram_budget"] + + +def test_fdtd_independent_emitters_scale_work_not_serial_memory(): + kwargs = dict(size_um=(10, 10, 1), cell_nm=10, duration_fs=100) + one = fdtd_estimate(**kwargs) + ensemble = fdtd_estimate(**kwargs, independent_source_runs=300) + assert ensemble["cell_updates_all_runs"] == 300 * one["cell_updates_all_runs"] + assert ensemble["raw_grid_bytes"] == one["raw_grid_bytes"] + double = fdtd_estimate(**kwargs, bytes_per_real=8) + assert double["raw_grid_bytes"] == 2 * one["raw_grid_bytes"] + assert "gpu_ideal_roofline_seconds_lower_bound" not in one + + +def test_fdtd_memory_fit_can_fail_and_cfl_uses_fastest_index(): + result = fdtd_estimate(size_um=(100, 100, 1), cell_nm=10, duration_fs=100, + hardware=Hardware(100, 20, 1000, 100, 16, 8)) + assert not result["fits_cpu_ram_budget"] + assert not result["fits_gpu_vram_budget"] + slower = fdtd_estimate(size_um=(1, 1, 1), cell_nm=10, duration_fs=100, minimum_index=2) + faster = fdtd_estimate(size_um=(1, 1, 1), cell_nm=10, duration_fs=100, minimum_index=1) + assert slower["courant_step_fs"] == 2 * faster["courant_step_fs"] + assert slower["time_steps_per_run"] <= faster["time_steps_per_run"] + + +@pytest.mark.parametrize("kwargs", [ + {"cell_nm": 0}, {"duration_fs": -1}, {"minimum_index": float("nan")}, + {"pml_cells": 0}, {"bytes_per_real": 3}, {"allocated_real_arrays_per_cell": 5}, + {"independent_source_runs": 1.5}, {"size_um": (1, 2)}, {"size_um": (1, 0, 1)}, +]) +def test_invalid_fdtd_assumptions(kwargs): + args = dict(size_um=(1, 1, 1), cell_nm=10, duration_fs=100) + args.update(kwargs) + with pytest.raises(ValueError): + fdtd_estimate(**args) + + +@pytest.mark.parametrize("value", [0, -1, 1.5, True]) +def test_invalid_planar_counts(value): + with pytest.raises(ValueError): + planar_estimate(wavelengths=value, finite_layers=2, emitter_depths=1, quadrature_evaluations=100) + + +def test_invalid_hardware_and_population(): + with pytest.raises(ValueError): + Hardware(0, 20, 1000, 100, 16, 8) + with pytest.raises(ValueError): + emitter_population(area_um2=1, depth_um=0, density_per_um3=1000) diff --git a/tests/test_spectral_eqe.py b/tests/test_spectral_eqe.py new file mode 100644 index 0000000..1db80c3 --- /dev/null +++ b/tests/test_spectral_eqe.py @@ -0,0 +1,359 @@ +import hashlib +import json +from unittest.mock import patch +import io + +import numpy as np +import pytest + +from spectral_eqe import ( + DispersiveStack, Layer, OpticalDataLake, OpticalRateBudget, + dipole_ensemble_eqe, random_dipoles, +) + + +TABLE = b"""REFERENCES: Synthetic test data, not a measured material +COMMENTS: Unit test only +DATA: + - type: tabulated nk + data: | + 0.4 1.5 0.01 + 0.6 1.7 0.03 +""" + + +def import_table(lake, name="film", raw=TABLE): + return lake.import_rii(name, raw, source="synthetic fixture", revision="test-v1") + + +def constant(lake, name, n, k=0): + raw = f"REFERENCES: Synthetic test\nDATA:\n - type: tabulated nk\n data: |\n 0.4 {n} {k}\n 0.6 {n} {k}\n".encode() + return import_table(lake, name, raw) + + +def test_catalog_csv_interpolation_and_integrity(tmp_path): + source = tmp_path / "input.csv" + source.write_text("wavelength_nm,n,k\n400,1.5,0.01\n600,1.7,0.03\n") + lake = OpticalDataLake(tmp_path / "lake") + material = lake.import_csv("film", source, source="measurement.csv", license="test-only", + citation="Synthetic", conditions="test") + np.testing.assert_allclose(material.index([400, 500, 600]), [1.5 + .01j, 1.6 + .02j, 1.7 + .03j]) + assert lake.list_datasets() == ["film"] + assert lake.load("film").metadata["conditions"] == "test" + digest = hashlib.sha256(source.read_bytes()).hexdigest() + (lake.root / f"{digest}.source").write_bytes(b"corrupt") + with pytest.raises(ValueError, match="Raw material checksum"): + lake.load("film") + + +def test_catalog_record_checksum_and_immutable_ids(tmp_path): + lake = OpticalDataLake(tmp_path) + import_table(lake) + import_table(lake) + with pytest.raises(FileExistsError): + constant(lake, "film", 2) + record = json.loads((tmp_path / "film.json").read_text()) + record["payload"]["metadata"]["citation"] = "tampered" + (tmp_path / "film.json").write_text(json.dumps(record)) + with pytest.raises(ValueError, match="Catalog checksum"): + lake.load("film") + + +@pytest.mark.parametrize("kind,pole", [("formula 1", ".1"), ("formula 2", ".01")]) +def test_rii_native_tables_and_formula(tmp_path, kind, pole): + lake = OpticalDataLake(tmp_path) + formula = f"REFERENCES: Synthetic formula\nDATA:\n - type: {kind}\n wavelength_range: 0.4 0.6\n coefficients: 0 1 {pole}\n".encode() + material = lake.import_rii("formula", formula, source="synthetic", revision="v1", assume_lossless=True) + expected = np.sqrt(1 + .5**2 / (.5**2 - .1**2)) + assert material.index([500])[0] == pytest.approx(expected + 0j) + assert material.metadata["license"] == "CC0-1.0" + assert material.assume_lossless + nk = import_table(lake) + assert nk.index([500])[0] == pytest.approx(1.6 + .02j) + + +def test_separate_n_and_k_tables_use_intersection(tmp_path): + raw = b"""REFERENCES: Synthetic separate grids +DATA: + - type: tabulated n + data: | + 0.4 1.5 + 0.6 1.7 + - type: tabulated k + data: | + 0.45 0.01 + 0.55 0.03 +""" + material = import_table(OpticalDataLake(tmp_path), raw=raw) + assert material.index([500])[0] == pytest.approx(1.6 + .02j) + with pytest.raises(ValueError, match="domain"): + material.index([400]) + + +def test_missing_k_requires_explicit_lossless_assumption(tmp_path): + raw = b"REFERENCES: Synthetic\nDATA:\n - type: tabulated n\n data: |\n .4 1.5\n .6 1.5\n" + material = import_table(OpticalDataLake(tmp_path), raw=raw) + with pytest.raises(ValueError, match="Extinction coefficient"): + material.index([500]) + + +def test_fetch_is_pinned_and_offline(tmp_path): + lake = OpticalDataLake(tmp_path) + commit = "a" * 40 + path = "database/data/main/Ag/nk/Johnson.yml" + with patch("spectral_eqe.urlopen", return_value=io.BytesIO(TABLE)) as remote: + material = lake.fetch_rii("ag", path, commit) + assert remote.call_args.args[0] == f"https://raw.githubusercontent.com/polyanskiy/refractiveindex.info-database/{commit}/{path}" + assert material.metadata["revision"] == commit + with patch("spectral_eqe.urlopen", side_effect=AssertionError("must be offline")): + assert lake.fetch_rii("ag", path, commit).index([500])[0] == pytest.approx(1.6 + .02j) + with pytest.raises(FileExistsError): + lake.fetch_rii("ag", path, "b" * 40) + + +@pytest.mark.parametrize("path,commit", [ + ("database/data/main/Ag/nk/Johnson.yml", "main"), + ("database/data/../private.yml", "a" * 40), + ("https://example.com/file.yml", "a" * 40), + ("database/data/a.yml?x=1", "a" * 40), + ("database\\data\\a.yml", "a" * 40), +]) +def test_fetch_rejects_unpinned_and_noncanonical_paths(tmp_path, path, commit): + with patch("spectral_eqe.urlopen", side_effect=AssertionError("no network")): + with pytest.raises(ValueError): + OpticalDataLake(tmp_path).fetch_rii("sample", path, commit) + + +def test_fetch_network_errors_are_not_hidden(tmp_path): + with patch("spectral_eqe.urlopen", side_effect=OSError("offline")): + with pytest.raises(OSError, match="offline"): + OpticalDataLake(tmp_path).fetch_rii("sample", "database/data/a.yml", "a" * 40) + + +@pytest.mark.parametrize("size", [4_000_000, 4_000_001]) +def test_fetch_dataset_size_boundary(tmp_path, size): + raw = TABLE + b"#" + b" " * (size - len(TABLE) - 1) + lake = OpticalDataLake(tmp_path) + with patch("spectral_eqe.urlopen", return_value=io.BytesIO(raw)): + if size > 4_000_000: + with pytest.raises(ValueError, match="4 MB"): + lake.fetch_rii("sample", "database/data/a.yml", "a" * 40) + assert lake.list_datasets() == [] + else: + material = lake.fetch_rii("sample", "database/data/a.yml", "a" * 40) + assert material.index([500])[0] == pytest.approx(1.6 + .02j) + + +def test_catalog_rejects_path_ids_and_missing_csv_provenance(tmp_path): + lake = OpticalDataLake(tmp_path / "lake") + with pytest.raises(ValueError, match="identifier"): + lake.load("../outside") + source = tmp_path / "data.csv" + source.write_text("wavelength_nm,n,k\n400,1,0\n600,1,0\n") + with pytest.raises(ValueError, match="license"): + lake.import_csv("sample", source, source="source", license="", citation="Test", conditions="Test") + assert lake.list_datasets() == [] + source.write_text("wavelength_um,n,k\n.4,1,0\n.6,1,0\n") + with pytest.raises(ValueError, match="columns"): + lake.import_csv("sample", source, source="source", license="test", citation="Test", conditions="Test") + + +def test_formula_domain_poles_and_malformed_coefficients_are_rejected(tmp_path): + lake = OpticalDataLake(tmp_path) + raw = b"REFERENCES: Test\nDATA:\n - type: formula 1\n wavelength_range: 0.4 0.6\n coefficients: 0 1 0.5\n" + material = lake.import_rii("pole", raw, source="fixture", revision="1", assume_lossless=True) + with pytest.raises(ValueError, match="pole"): + material.index([500]) + with pytest.raises(ValueError, match="domain"): + material.index([700]) + with pytest.raises(ValueError): + lake.import_rii("malformed", raw.replace(b"0 1 0.5", b"0 1 nonsense"), source="fixture", revision="1") + + +def test_cached_lossless_assumption_cannot_change(tmp_path): + lake = OpticalDataLake(tmp_path) + with patch("spectral_eqe.urlopen", return_value=io.BytesIO(TABLE)): + lake.fetch_rii("sample", "database/data/a.yml", "a" * 40) + with pytest.raises(FileExistsError, match="assumption"): + lake.fetch_rii("sample", "database/data/a.yml", "a" * 40, assume_lossless=True) + + +@pytest.mark.parametrize("raw", [ + b"REFERENCES: Test\nDATA: []", + b"REFERENCES: Test\nDATA:\n - type: formula 9\n", + b"DATA:\n - type: tabulated nk\n data: |\n .4 1.5 0\n .6 1.5 0\n", + b"REFERENCES: Test\nDATA:\n - type: tabulated nk\n data: |\n .4 1.5 -1\n .6 1.5 0\n", + b"REFERENCES: Test\nDATA:\n - type: tabulated nk\n data: |\n .6 1.5 0\n .4 1.5 0\n", +]) +def test_invalid_imports_are_explicit(tmp_path, raw): + with pytest.raises(ValueError): + import_table(OpticalDataLake(tmp_path), raw=raw) + + +@pytest.mark.parametrize("grid", [[399], [601], [500, 400], [500, 500], [], [np.nan], [[500]]]) +def test_material_rejects_invalid_or_out_of_domain_wavelengths(tmp_path, grid): + material = import_table(OpticalDataLake(tmp_path)) + with pytest.raises(ValueError): + material.index(grid) + + +def make_stack(lake, thickness=80): + constant(lake, "air", 1) + constant(lake, "emitter", 1.7) + import_table(lake) + return DispersiveStack((Layer("air", None), Layer("emitter", 40), + Layer("film", thickness), Layer("air", None)), 1, 20) + + +def test_stack_dispersion_and_thickness(tmp_path): + lake = OpticalDataLake(tmp_path) + stack = make_stack(lake) + data = stack.evaluate(lake, [400, 500, 600]) + assert data["indices"].shape == (3, 4) + np.testing.assert_allclose(data["indices"][:, 2], [1.5 + .01j, 1.6 + .02j, 1.7 + .03j]) + np.testing.assert_equal(data["thicknesses_nm"], [np.inf, 40, 80, np.inf]) + assert data["emitter_depth_nm"] == 20 + assert data["materials"][2]["dataset_id"] == "film" + + +def test_plane_wave_thickness_spectrum_is_not_eqe(tmp_path): + lake = OpticalDataLake(tmp_path) + thin = make_stack(lake, 10).plane_wave_spectrum(lake, [400, 500, 600]) + thick = make_stack(lake, 100).plane_wave_spectrum(lake, [400, 500, 600]) + np.testing.assert_allclose(thin["R"] + thin["T"] + thin["A"], 1, atol=1e-12) + assert thin["quantity"] == "plane_wave_power_not_EQE" + assert np.max(np.abs(thin["T"] - thick["T"])) > .01 + assert np.all(thin["A"] >= 0) + assert "EQE" not in thin + + +@pytest.mark.parametrize("layers,emitter,depth", [ + ((Layer("a", None), Layer("b", None)), 1, 1), + ((Layer("a", 1), Layer("b", 10), Layer("c", None)), 1, 1), + ((Layer("a", None), Layer("b", -1), Layer("c", None)), 1, 1), + ((Layer("a", None), Layer("b", 10), Layer("c", None)), 0, 1), + ((Layer("a", None), Layer("b", 10), Layer("c", None)), 1, 0), + ((Layer("a", None), Layer("b", 10), Layer("c", None)), 1, 10), + ((Layer("a", None), Layer("b", 10), Layer("c", None)), 1, np.nan), +]) +def test_invalid_stack_geometry(layers, emitter, depth): + with pytest.raises(ValueError): + DispersiveStack(layers, emitter, depth) + + +def test_absorbing_host_and_halfspace_are_rejected(tmp_path): + lake = OpticalDataLake(tmp_path) + constant(lake, "air", 1) + constant(lake, "lossy", 1.7, .1) + for layers in [ + (Layer("air", None), Layer("lossy", 10), Layer("air", None)), + (Layer("lossy", None), Layer("air", 10), Layer("air", None)), + ]: + with pytest.raises(ValueError, match="lossless"): + DispersiveStack(layers, 1, 5).evaluate(lake, [500]) + + +def test_monochromatic_rate_equation_and_yield_boundaries(): + budget = OpticalRateBudget([400, 500], [1, 2], [.5, .4], "P_host") + kwargs = dict(charge_balance=.8, exciton_yield=.25) + np.testing.assert_allclose(budget.monochromatic_eqe(intrinsic_yield=.5, **kwargs), [.05, .2 * .2 / 1.5]) + np.testing.assert_array_equal(budget.monochromatic_eqe(intrinsic_yield=0, **kwargs), [0, 0]) + np.testing.assert_allclose(budget.monochromatic_eqe(intrinsic_yield=1, **kwargs), [.1, .04]) + np.testing.assert_array_equal( + budget.monochromatic_eqe(intrinsic_yield=1, charge_balance=0, exciton_yield=1), [0, 0]) + + +def test_broadband_uses_global_rate_denominator(): + budget = OpticalRateBudget([400, 500, 600], [1, 2, 3], [.5, 1, 1.5], "P_host") + result = budget.broadband_eqe([.005, .005, .005], intrinsic_yield=.5, charge_balance=.8, exciton_yield=.25) + assert result["rate_denominator"] == pytest.approx(1.5) + np.testing.assert_allclose(result["eqe_density_per_nm"], .2 * .5 * .005 * np.array([.5, 1, 1.5]) / 1.5) + assert result["integrated_eqe"] == pytest.approx(.2 * .5 / 1.5) + zero = budget.broadband_eqe([.005] * 3, intrinsic_yield=0, charge_balance=1, exciton_yield=1) + assert zero["integrated_eqe"] == 0 + np.testing.assert_array_equal(zero["eqe_density_per_nm"], [0, 0, 0]) + perfect = budget.broadband_eqe([.005] * 3, intrinsic_yield=1, charge_balance=1, exciton_yield=1) + assert perfect["integrated_eqe"] == pytest.approx(.5) + + +def test_orientation_averages_powers_before_ratios(): + vertical = OpticalRateBudget([500], [4], [1], "P_host") + horizontal = OpticalRateBudget([500], [1], [.5], "P_host") + average = OpticalRateBudget.orientation_average(vertical, horizontal, 1 / 3) + assert average.total[0] == pytest.approx(2) + assert average.escape[0] == pytest.approx(2 / 3) + assert average.monochromatic_eqe(intrinsic_yield=1, charge_balance=1, exciton_yield=1)[0] == pytest.approx(1 / 3) + compare = dipole_ensemble_eqe(vertical, horizontal, [[0, 0, 1], [1, 0, 0]], + intrinsic_yield=1, charge_balance=1, exciton_yield=1) + assert compare["independent_fixed_dipoles_eqe"][0] == pytest.approx((.25 + .5) / 2) + assert compare["rapid_reorientation_eqe"][0] == pytest.approx(.75 / 2.5) + assert compare["mean_vertical_fraction"] == .5 + + +def test_random_orientations_are_reproducible_and_isotropic(): + a = random_dipoles(10000, seed=42) + np.testing.assert_array_equal(a, random_dipoles(10000, seed=42)) + np.testing.assert_allclose(np.linalg.norm(a, axis=1), 1, atol=1e-14) + np.testing.assert_allclose(np.mean(a**2, axis=0), 1 / 3, atol=.01) + assert not np.array_equal(a, random_dipoles(10000, seed=43)) + + +@pytest.mark.parametrize("total,escape", [([-1], [0]), ([1], [-1]), ([1], [2]), ([np.nan], [0]), ([1, 2], [.5])]) +def test_invalid_optical_budgets(total, escape): + with pytest.raises(ValueError): + OpticalRateBudget([500], total, escape, "P_host") + + +@pytest.mark.parametrize("yield_value", [-1, 1.1, np.nan, np.inf]) +def test_invalid_eqe_factors(yield_value): + budget = OpticalRateBudget([500], [1], [.5], "P_host") + for parameter in ["intrinsic_yield", "charge_balance", "exciton_yield"]: + kwargs = dict(intrinsic_yield=1, charge_balance=1, exciton_yield=1) + kwargs[parameter] = yield_value + with pytest.raises(ValueError): + budget.monochromatic_eqe(**kwargs) + + +def test_dark_emitter_and_invalid_spectrum_fail_explicitly(): + budget = OpticalRateBudget([400, 600], [0, 0], [0, 0], "P_host") + with pytest.raises(ValueError, match="Zero total"): + budget.monochromatic_eqe(intrinsic_yield=1, charge_balance=1, exciton_yield=1) + for spectrum in [[1, 1], [-.005, .015], [np.nan, .005], [.005]]: + with pytest.raises(ValueError): + budget.broadband_eqe(spectrum, intrinsic_yield=.5, charge_balance=1, exciton_yield=1) + with pytest.raises(ValueError, match="Zero integrated"): + budget.broadband_eqe([.005, .005], intrinsic_yield=1, charge_balance=1, exciton_yield=1) + + +def test_spectral_overflow_is_not_reported_as_zero_efficiency(): + # A narrow normalized peak exposes intermediate spectrum-times-rate overflow. + budget = OpticalRateBudget([500, 500.0000001], [1.7e308] * 2, [1.7e308] * 2, "P_host") + density = 1 / (500.0000001 - 500) + with pytest.raises(ValueError, match="overflowed"): + budget.broadband_eqe([density] * 2, intrinsic_yield=1, charge_balance=1, exciton_yield=1) + + +def test_budget_copies_inputs_and_requires_matching_references(): + total = np.array([1.]) + budget = OpticalRateBudget([500], total, [.5], "host A") + total[0] = 2 + assert budget.total[0] == 1 + with pytest.raises(ValueError): + budget.total[0] = 2 + for other in [OpticalRateBudget([500], [1], [.5], "host B"), + OpticalRateBudget([600], [1], [.5], "host A")]: + with pytest.raises(ValueError): + OpticalRateBudget.orientation_average(budget, other, .5) + + +@pytest.mark.parametrize("directions", [[], [[0, 0, 0]], [[np.nan, 0, 1]], [[1, 0]]]) +def test_invalid_ensemble_directions(directions): + budget = OpticalRateBudget([500], [1], [.5], "P_host") + with pytest.raises(ValueError): + dipole_ensemble_eqe(budget, budget, directions, intrinsic_yield=1, charge_balance=1, exciton_yield=1) + + +@pytest.mark.parametrize("count", [0, -1, True, 1.5]) +def test_random_orientation_count_validation(count): + with pytest.raises(ValueError): + random_dipoles(count, seed=42) From 6c7738c9f7c6370b3eb87b2d3cac52f2c5de53b2 Mon Sep 17 00:00:00 2001 From: "Guy Francis Mongelli, MS, EIT" Date: Fri, 16 Oct 2026 18:48:52 -0400 Subject: [PATCH 6/7] docs: preserve full mathematical derivations with rendered equations Consolidate dipole source work, dispersion, spectral EQE, orientation populations and resource-planning equations with rationales, assumptions, references and explicit implementation boundaries. Link the record from companion guides and include it in the source distribution. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- MANIFEST.in | 1 + README.md | 4 + docs/DIPOLE_FORMULATION.md | 3 + docs/ENSEMBLE_RESOURCE_ESTIMATES.md | 3 + docs/MATHEMATICAL_DERIVATIONS.md | 853 ++++++++++++++++++++++++++++ docs/SPECTRAL_EQE_AND_DATA.md | 3 + 6 files changed, 867 insertions(+) create mode 100644 docs/MATHEMATICAL_DERIVATIONS.md diff --git a/MANIFEST.in b/MANIFEST.in index 4576631..dda94aa 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,3 +1,4 @@ include docs/DIPOLE_FORMULATION.md include docs/SPECTRAL_EQE_AND_DATA.md include docs/ENSEMBLE_RESOURCE_ESTIMATES.md +include docs/MATHEMATICAL_DERIVATIONS.md diff --git a/README.md b/README.md index 4a36da3..147b096 100644 --- a/README.md +++ b/README.md @@ -127,6 +127,10 @@ license presentation, final screenshots, seller verification and publication. ### Dispersive data, spectral EQE and ensemble planning +The [consolidated mathematical derivations](docs/MATHEMATICAL_DERIVATIONS.md) +render full LaTeX equations on GitHub, with assumptions, step-by-step rationales, +benchmark limits, resource models and implementation-status boundaries. + See [native optical data and spectral EQE derivation](docs/SPECTRAL_EQE_AND_DATA.md) and [orientation/resource planning](docs/ENSEMBLE_RESOURCE_ESTIMATES.md). The new APIs provide pinned native refractiveindex.info ingestion, local diff --git a/docs/DIPOLE_FORMULATION.md b/docs/DIPOLE_FORMULATION.md index 8f4a5a5..58a4e4b 100644 --- a/docs/DIPOLE_FORMULATION.md +++ b/docs/DIPOLE_FORMULATION.md @@ -1,5 +1,8 @@ # Dipole reference kernels: conventions and derivation +For rendered LaTeX equations and the connected EQE/ensemble/resource derivations, +see the [consolidated mathematical record](MATHEMATICAL_DERIVATIONS.md). + Research milestone, 2026-10-16. Implemented in `dipole_reference.py`. These are homogeneous and **one-interface** reference kernels, not a finite emitting-layer/multilayer OLED solver, material database, or sale-ready product. diff --git a/docs/ENSEMBLE_RESOURCE_ESTIMATES.md b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md index 3e09b08..91430f5 100644 --- a/docs/ENSEMBLE_RESOURCE_ESTIMATES.md +++ b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md @@ -1,5 +1,8 @@ # Orientation ensembles and CPU/GPU/RAM planning +See the [consolidated mathematical record](MATHEMATICAL_DERIVATIONS.md) +for rendered LaTeX equations, assumptions and connected optical/EQE derivations. + Prepared 2026-10-16. APIs: `spectral_eqe.py` and `simulation_resources.py`. Estimates are explicit workload models, NOT a GPU implementation, hardware performance guarantee, measured peak memory, or full-pixel field simulation. diff --git a/docs/MATHEMATICAL_DERIVATIONS.md b/docs/MATHEMATICAL_DERIVATIONS.md new file mode 100644 index 0000000..42d3bfe --- /dev/null +++ b/docs/MATHEMATICAL_DERIVATIONS.md @@ -0,0 +1,853 @@ +# OLED optics, spectral EQE and ensemble-resource derivations + +**Research record: 2026-10-16.** This document uses GitHub-supported LaTeX +math, rendered directly in the repository. Open the rendered file on GitHub +or in a Markdown preview with math support, rather than a plain-text viewer. +Equations are restated with project-specific explanations; no third-party +chapter text or figures are reproduced. + +## 1. Scope and reading map + +This is the consolidated mathematical record for the implemented foundations, +their rationale, numerical transformations and planning estimates. + +| Subject | Status | Implementation / companion | +|---|---|---| +| Homogeneous angular emission and one-interface source work | Implemented, benchmarked reference kernels | [dipole_reference.py](../dipole_reference.py), [formulation](DIPOLE_FORMULATION.md) | +| Native optical data, dispersion and layer geometry | Implemented | [spectral_eqe.py](../spectral_eqe.py), [data guide](SPECTRAL_EQE_AND_DATA.md) | +| Dispersive incident-plane-wave R/T/A | Implemented diagnostic, not emitted flux | [spectral_eqe.py](../spectral_eqe.py) | +| Monochromatic/broadband EQE from supplied physical rate budgets | Implemented composition, not automatic stack-derived EQE | [spectral_eqe.py](../spectral_eqe.py) | +| Fixed versus reorienting dipole ensembles | Implemented for monochromatic transition arrays | [spectral_eqe.py](../spectral_eqe.py) | +| Planar workload, finite-pixel FDTD memory/work and optional roofline | Planning models only | [simulation_resources.py](../simulation_resources.py), [resource guide](ENSEMBLE_RESOURCE_ESTIMATES.md) | +| Finite emitting-layer multilayer source/escape calculation | Not implemented | [solver acceptance specification](OLED_SOLVER_SPEC.md) | +| Spatially integrated device EQE, finite-pixel FDTD, GPU backend | Not implemented | Future work; formulas below identify requirements, not delivered solvers | + +In particular, source work, external illumination transmission, collected +emitted power and EQE are different observables. A mathematically valid +composition formula does not validate its supplied material data or rates. + +## 2. Conventions, units and power normalization + +Assume linear, isotropic, nonmagnetic media and harmonic peak amplitudes with +time dependence + +$$ +\mathbf E(\mathbf r,t)=\operatorname{Re} + \left[\mathbf E(\mathbf r)e^{-i\omega t}\right], +\qquad +k_0=\frac{2\pi}{\lambda_0}=\frac{\omega}{c}. +$$ + +The relative index and permittivity are + +$$ +N_j(\lambda_0)=n_j(\lambda_0)+i\kappa_j(\lambda_0), +\qquad \epsilon_j=N_j^2. +$$ + +The symbol $\kappa$ denotes extinction coefficient, not a wavevector. +Passive data have $\kappa_j\geq0$ and, in the supported index convention, +$n_j\geq0$. The emitting reference host has real $n_h>0$ and +$k_h=n_hk_0$. Fields propagate as $e^{ik_z z}$ in their outgoing coordinate; +choose longitudinal roots that propagate outward or decay rather than grow. +The API uses vacuum wavelengths and geometric lengths in nm. Ratios and +phases must use consistent length units; SI units are required in power +formulas with $c$, $\epsilon_0$ and a dipole moment in C m. + +For a prescribed electric dipole $\mathbf p$, the homogeneous-host power is + +$$ +P_h=\frac{n_h\omega^4|\mathbf p|^2}{12\pi\epsilon_0c^3}. +\tag{1} +$$ + +This is a macroscopic lossless-host normalization for a **peak** amplitude, +not an RMS amplitude. No microscopic local-field correction is included. +Using the same dipole and host reference throughout is essential when +combining rates from different orientations or collection channels. + +The field radiated back onto the dipole changes source work. With the +outgoing electric dyadic Green tensor defined by +$\mathbf E_{\rm ref}=\mu_0\omega^2\mathbf G_{\rm ref}\mathbf p$, + +$$ +P_{\rm ref}=\frac{\omega}{2} + \operatorname{Im}\!\left[\mathbf p^*\cdot\mathbf E_{\rm ref}\right], +\qquad +F_{\rm EM}=1+\frac{6\pi}{k_h} + \operatorname{Im}\!\left[ + \hat{\mathbf e}^{\,*}\cdot\mathbf G_{\rm ref}\hat{\mathbf e}\right], +\tag{2} +$$ + +where $\hat{\mathbf e}=\mathbf p/|\mathbf p|$. The coefficient follows by +dividing the source-work expression by (1), using +$\mu_0\epsilon_0c^2=1$. The Green tensor has units of inverse length. +For absorbing surroundings, source work includes dissipated electromagnetic +power as well as outgoing radiation. It excludes intrinsic emitter decay +that is not represented by the electromagnetic environment. + +## 3. Homogeneous angular pattern and polarization + +For observation direction $\hat{\mathbf r}$, the transverse dipole field is +proportional to +$\hat{\mathbf r}\times(\hat{\mathbf e}\times\hat{\mathbf r})$. +Its squared magnitude is $1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2$. +Normalizing its solid-angle integral gives + +$$ +\frac{1}{P_h}\frac{dP}{d\Omega} + =\frac{3}{8\pi} + \left(1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2\right), +\qquad +d\Omega=\sin\theta\,d\theta\,d\phi. +\tag{3} +$$ + +Indeed, $\int d\Omega=4\pi$ and +$\int|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2d\Omega=4\pi/3$, +so the normalized integral is one. Reflection symmetry gives one-half +in either homogeneous hemisphere. + +Project the dipole onto the orthogonal far-field polarization axes before +squaring: + +$$ +\frac{dP_s}{P_h\,d\Omega}=\frac{3}{8\pi} + |\hat{\mathbf e}\cdot\hat{\boldsymbol\phi}|^2, +\qquad +\frac{dP_p}{P_h\,d\Omega}=\frac{3}{8\pi} + |\hat{\mathbf e}\cdot\hat{\boldsymbol\theta}|^2. +\tag{4} +$$ + +| Dipole distribution | $s$ factor, excluding $3/(8\pi)$ | $p$ factor, excluding $3/(8\pi)$ | +|---|---|---| +| Vertical, $z$ | $0$ | $\sin^2\theta$ | +| Horizontal, $x$ | $\sin^2\phi$ | $\cos^2\theta\cos^2\phi$ | +| Horizontal, $y$ | $\cos^2\phi$ | $\cos^2\theta\sin^2\phi$ | +| Uniform random in-plane | $1/2$ | $\cos^2\theta/2$ | +| Uniform isotropic | $1/3$ | $1/3$ | + +These averages are averages of **powers**. Averaging unrelated complex +dipole amplitudes first would incorrectly cancel incoherent emission. + +## 4. One-interface source-work kernels + +Place the interface at $z=0$ and the dipole at $z=d>0$ in medium 1, +with $\epsilon_1=n_h^2$. Medium 2 occupies $z<0$. +For each in-plane wavevector magnitude, define + +$$ +q=\frac{k_\parallel}{k_h},\qquad +w=\sqrt{1-q^2},\qquad +v=\sqrt{\frac{\epsilon_2}{\epsilon_1}-q^2}, +\qquad a=2k_hd=\frac{4\pi n_hd}{\lambda_0}. +\tag{5} +$$ + +On the real $q$ axis, $w$ is positive real below the host light line +and positive imaginary above it. The substrate root uses the corresponding +outgoing/decaying prescription. Electric-field Fresnel reflection amplitudes +are + +$$ +r_s=\frac{w-v}{w+v},\qquad +r_p=\frac{\epsilon_2w-\epsilon_1v} + {\epsilon_2w+\epsilon_1v}. +\tag{6} +$$ + +Identical media give $r_s=r_p=0$. The explicit perfect-electric-conductor +boundary uses $r_s=-1$, $r_p=1$; it is an ideal boundary, not a finite +tabulated metal index. + +The angular-spectrum representation of the reflected Green tensor, +followed by azimuthal integration and division by (1), gives + +$$ +F_\perp=1+\frac32\operatorname{Re}\int_0^\infty + \frac{q^3}{w}r_p e^{iaw}\,dq, +\tag{7} +$$ + +$$ +F_\parallel=1+\frac34\operatorname{Re}\int_0^\infty + \frac{q}{w}\left(r_s-w^2r_p\right)e^{iaw}\,dq. +\tag{8} +$$ + +Here $\perp$ means dipole normal to the interface, not the $s$ polarization. +The phase is a round trip from source to interface and back. +The **minus sign** in $r_s-w^2r_p$ comes from the reflected $p$-polarization +dyad: incident and reflected polarization vectors have different longitudinal +directions. An intensity average cannot replace this signed field product. +The ideal-mirror tests below independently constrain this sign. + +The $q<1$ / $q>1$ split classifies propagation **in the host**. +A host-evanescent component can propagate in a higher-index substrate. +Thus neither the evanescent correction nor $F_{\rm EM}$ alone is an +absorption fraction or an external collection efficiency. + +### 4.1 Removing the host light-line singularity + +At $q=1$, the $1/w$ factor is integrable but unsuitable for direct point +sampling. In the propagating sector use $u=\sqrt{1-q^2}$, so +$q\,dq=-u\,du$. Reversing the integration limits yields + +$$ +\Delta F_{\perp,\rm prop}=\frac32\int_0^1 + (1-u^2)\operatorname{Re}[r_pe^{iau}]\,du, +\tag{9} +$$ + +$$ +\Delta F_{\parallel,\rm prop}=\frac34\int_0^1 + \operatorname{Re}[(r_s-u^2r_p)e^{iau}]\,du. +\tag{10} +$$ + +For the evanescent sector set $t=\sqrt{q^2-1}$, then $w=it$, +$q\,dq=t\,dt$ and $e^{iaw}=e^{-at}$. Finally let $x=at$: + +$$ +\Delta F_{\perp,\rm evan}=\frac{3}{2a}\int_0^\infty + \left[1+\left(\frac{x}{a}\right)^2\right] + e^{-x}\operatorname{Im}r_p\,dx, +\tag{11} +$$ + +$$ +\Delta F_{\parallel,\rm evan}=\frac{3}{4a}\int_0^\infty + e^{-x}\operatorname{Im}\left[ + r_s+\left(\frac{x}{a}\right)^2r_p\right]dx. +\tag{12} +$$ + +The coefficients in each transformed integral are evaluated at that integral's +$q$. The Jacobians cancel the host light-line singularity; scaling by $a$ +also makes the dominant decay envelope $e^{-x}$ independent of separation. +The full answer is $1+\Delta F_{\rm prop}+\Delta F_{\rm evan}$. + +### 4.2 Branches, resonances and error gates + +Write $\rho=\epsilon_2/\epsilon_1$. For a lossless dielectric, its light line +appears at + +$$ +u_b=\sqrt{1-\rho}\quad(0<\rho<1), +\qquad x_b=a\sqrt{\rho-1}\quad(\rho>1). +\tag{13} +$$ + +Splitting there helps resolve changes in the substrate root. Weak loss moves +these branch locations off the real axis; the integration intervals are split +around their projected centers and widths. + +A $p$ surface-mode denominator in (6) gives +$q_{\rm sp}^2=\rho/(\rho+1)$ and hence + +$$ +t_{\rm sp}=\sqrt{-\frac{1}{\rho+1}},\qquad x_{\rm sp}=at_{\rm sp}. +\tag{14} +$$ + +For supported lossy-metal cases, intervals around +$\operatorname{Re}x_{\rm sp}$ span geometrically increasing multiples of +$|\operatorname{Im}x_{\rm sp}|$. This resolves narrow resonances and wings +better than a single large adaptive interval. + +The implementation integrates to infinity, using an initial tail split +$x=32$ and a refinement split $x=64$. Their quadrature tolerance factors +are $0.05$ and $0.0125$ of the requested tolerances, respectively. +Per-interval absolute budgets share the scaled absolute tolerance. +For each orientation it checks the summed quadrature error estimate and + +$$ +\delta_{\rm refine}= + |\Delta F_{\rm prop}^{(2)}-\Delta F_{\rm prop}^{(1)}| + +|\Delta F_{\rm evan}^{(2)}-\Delta F_{\rm evan}^{(1)}| + \leq \max(\mathrm{atol},\mathrm{rtol}|F|). +\tag{15} +$$ + +Failures, nonfinite results and negative total rates raise errors; rates are +not clipped into plausible values. The reported absolute signed tail integral +is a diagnostic, not a rigorous error bound. The exponential integrand is +evaluated as zero beyond $x=700$, outside supported floating-point exponential +resolution at the requested accuracy. + +The current kernel rejects $a<10^{-6}$, absorbing hosts and strictly lossless +nonpositive-permittivity substrates. Real-axis poles and ENZ limits need a +separate limiting-absorption prescription. Numerical support does not imply +that local continuum material constants remain valid at molecular distances. + +### 4.3 Independent analytic limits + +For an ideal mirror, image-dipole interference gives + +$$ +F_\perp=1+3\frac{\sin a-a\cos a}{a^3}, +\qquad +F_\parallel=1-\frac32\left( + \frac{\sin a}{a}+\frac{\cos a}{a^2}-\frac{\sin a}{a^3}\right). +\tag{16} +$$ + +Series expansion avoids cancellation near the surface: + +$$ +F_\perp=2-\frac{a^2}{10}+\frac{a^4}{280}+O(a^6), +\qquad +F_\parallel=\frac{a^2}{5}-\frac{3a^4}{280}+O(a^6). +\tag{17} +$$ + +The vertical image reinforces the source while the horizontal image cancels +it at vanishing separation. Both enhancements tend to one far from the +interface. These different limits detect sign, phase and prefactor mistakes. + +For a lossy half-space at small $a$, large $q$ dominates: +$r_s\to0$, $r_p\to(\epsilon_2-\epsilon_1)/(\epsilon_2+\epsilon_1)$. +Using $\int_0^\infty x^2e^{-x}dx=2$ in (11)-(12) yields + +$$ +\Delta F_\perp\sim\frac{3}{a^3} + \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}, +\qquad +\Delta F_\parallel\sim\frac{3}{2a^3} + \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}. +\tag{18} +$$ + +This explains the leading inverse-cube dependence and two-to-one orientation +ratio. It is an asymptotic check, not a substitute for converged integration +or an assertion of material validity at arbitrarily small distances. + +## 5. Dispersion, native data and layer thickness + +For tabulated data between $\lambda_i$ and $\lambda_{i+1}$, the current +interpolator applies independently to $n$ and $\kappa$: + +$$ +f(\lambda)=f_i+ + \frac{\lambda-\lambda_i}{\lambda_{i+1}-\lambda_i}(f_{i+1}-f_i), +\qquad f\in\{n,\kappa\}. +\tag{19} +$$ + +Native database wavelengths in micrometers are converted to nm: +$\lambda_{\rm nm}=1000\lambda_{\rm um}$. +No extrapolation is allowed. Separate $n$ and $\kappa$ tables must both cover +the query, so their usable domain is their intersection. Missing extinction +data require an explicit, recorded lossless assumption; missing data are not +evidence that absorption vanishes. + +The supported native Sellmeier records use $\lambda$ in **micrometers**: + +$$ +\begin{aligned} +\text{formula 1:}\quad +n^2(\lambda)=1+c_0+\sum_m + \frac{B_m\lambda^2}{\lambda^2-C_m^2},\\ +\text{formula 2:}\quad +n^2(\lambda)=1+c_0+\sum_m + \frac{B_m\lambda^2}{\lambda^2-C_m}. +\end{aligned} +\tag{20} +$$ + +The two pole conventions must not be interchanged. Evaluation must remain +inside the source range, avoid poles and produce finite positive $n^2$. +Unsupported formulas are rejected, not approximated by these forms. + +Each layer has its own $N_j(\lambda)$ and thickness $d_j$. For a plane-wave +component with conserved $k_\parallel$, + +$$ +k_{z,j}=\sqrt{N_j^2k_0^2-k_\parallel^2}, +\qquad \delta_j=k_{z,j}d_j. +\tag{21} +$$ + +Thickness changes phase and attenuation, while dispersion changes propagation +and Fresnel coefficients. They therefore affect multilayer interference jointly, +not as independent scalar correction factors. + +The implemented diagnostic calls the coherent transfer-matrix calculation +for both incident polarizations and averages their **powers**: + +$$ +R=\frac{R_s+R_p}{2},\qquad +T=\frac{T_s+T_p}{2},\qquad A=1-R-T. +\tag{22} +$$ + +For the enforced lossless exterior media, power transmission includes +the appropriate flux/admittance ratio, not just $|t|^2$. +This accounts for externally incident illumination. It is **not** the +collected power of a dipole inside the stack, whose source terms and +multiple-reflection feedback differ. + +The catalog preserves original source bytes, a pinned source revision, +license/citation/conditions and SHA-256 integrity checks. Checksums identify +changes, not physical accuracy or trusted authorship. The official +refractiveindex.info database is CC0; that does not license unrelated +databases or publication figures/text. See the [data guide](SPECTRAL_EQE_AND_DATA.md) +for the pinned Ag/SiO2 example and provenance details. + +## 6. From electromagnetic rates to EQE + +### 6.1 Monochromatic transition + +Define an intrinsic radiative rate in the homogeneous reference, +$\Gamma_{r,h}$, and a fixed intrinsic nonradiative rate $\Gamma_{nr}$: + +$$ +q_0=\frac{\Gamma_{r,h}}{\Gamma_{r,h}+\Gamma_{nr}}, +\qquad +\Gamma_{nr}=\Gamma_{r,h}\frac{1-q_0}{q_0} +\quad(0 Date: Sat, 17 Oct 2026 17:37:26 -0400 Subject: [PATCH 7/7] docs: compact equations and use operatorname macros Collapse multi-line display equations, convert plaintext equations in the companion docs to GitHub LaTeX, and standardize function/label macros on operatorname. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- docs/DIPOLE_FORMULATION.md | 57 +++--- docs/ENSEMBLE_RESOURCE_ESTIMATES.md | 15 +- docs/MATHEMATICAL_DERIVATIONS.md | 275 +++++++--------------------- docs/SPECTRAL_EQE_AND_DATA.md | 27 ++- 4 files changed, 127 insertions(+), 247 deletions(-) diff --git a/docs/DIPOLE_FORMULATION.md b/docs/DIPOLE_FORMULATION.md index 58a4e4b..e052110 100644 --- a/docs/DIPOLE_FORMULATION.md +++ b/docs/DIPOLE_FORMULATION.md @@ -19,13 +19,17 @@ ideal-boundary sentinel `"pec"` for a perfect electric conductor. Homogeneous macroscopic radiated power, for peak (not RMS) dipole amplitude: - P_host = n1 * omega^4 * |p|^2 / (12 pi epsilon0 c^3) +$$ +P_{\operatorname{host}}=\frac{n_1\omega^4|\mathbf p|^2}{12\pi\epsilon_0c^3} +$$ No microscopic local-field correction is included. The normalized angular density is dP/(P_host dOmega), with dOmega=sin(theta) dtheta dphi. For a unit dipole axis e and observation direction r_hat: - density = 3/(8 pi) * (1 - |e . r_hat|^2) +$$ +\frac{1}{P_{\operatorname{host}}}\frac{dP}{d\Omega}=\frac{3}{8\pi}\left(1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2\right) +$$ Project e onto e_phi (s) and e_theta (p) before squaring. In particular: @@ -48,17 +52,19 @@ v=sqrt(epsilon2/epsilon1-q^2). Select the outgoing/decaying branch: positive real root for propagation, positive imaginary root for evanescence. The p coefficient here uses the standard electric-field Fresnel convention: - r_s = (w-v)/(w+v) - r_p = (epsilon2*w-epsilon1*v)/(epsilon2*w+epsilon1*v) +$$ +r_s=\frac{w-v}{w+v},\qquad r_p=\frac{\epsilon_2w-\epsilon_1v}{\epsilon_2w+\epsilon_1v} +$$ For `"pec"`, r_s=-1 and r_p=+1. Identical media give both coefficients zero. The total electromagnetic source-work enhancement relative to P_host is: - F_vertical = 1 + (3/2) Re integral_0^infinity - [q^3/w * r_p * exp(2 i k1 d w)] dq - - F_parallel = 1 + (3/4) Re integral_0^infinity - [q/w * (r_s - w^2 r_p) * exp(2 i k1 d w)] dq +$$ +\begin{aligned} +F_\perp&=1+\frac32\operatorname{Re}\int_0^\infty\frac{q^3}{w}\,r_p\,e^{2ik_1dw}\,dq\\ +F_\parallel&=1+\frac34\operatorname{Re}\int_0^\infty\frac{q}{w}\,(r_s-w^2r_p)\,e^{2ik_1dw}\,dq +\end{aligned} +$$ The **minus** before w^2*r_p in the parallel expression is essential. It comes from the reflected p dyad, not from a polarization-intensity average. @@ -76,17 +82,21 @@ lossless substrate; "evanescent correction" does not mean "absorbed power". Write a=2 k1 d. In the q<1 region substitute u=sqrt(1-q^2): - Delta_vertical_prop = (3/2) integral_0^1 - (1-u^2) Re[r_p exp(i a u)] du - Delta_parallel_prop = (3/4) integral_0^1 - Re[(r_s-u^2*r_p) exp(i a u)] du +$$ +\begin{aligned} +\Delta F_{\perp,\operatorname{prop}}&=\frac32\int_0^1(1-u^2)\operatorname{Re}\left[r_pe^{iau}\right]du\\ +\Delta F_{\parallel,\operatorname{prop}}&=\frac34\int_0^1\operatorname{Re}\left[(r_s-u^2r_p)e^{iau}\right]du +\end{aligned} +$$ In q>1 substitute t=sqrt(q^2-1), then x=a*t: - Delta_vertical_evan = (3/2)/a integral_0^infinity - (1+(x/a)^2) exp(-x) Im[r_p] dx - Delta_parallel_evan = (3/4)/a integral_0^infinity - exp(-x) Im[r_s+(x/a)^2*r_p] dx +$$ +\begin{aligned} +\Delta F_{\perp,\operatorname{evan}}&=\frac{3}{2a}\int_0^\infty\left(1+\frac{x^2}{a^2}\right)e^{-x}\operatorname{Im}r_p\,dx\\ +\Delta F_{\parallel,\operatorname{evan}}&=\frac{3}{4a}\int_0^\infty e^{-x}\operatorname{Im}\left[r_s+\frac{x^2}{a^2}r_p\right]dx +\end{aligned} +$$ This removes the integrable q=1 singularity and scales the decaying tail. For a lossless dielectric, split the substrate light line explicitly: @@ -120,8 +130,12 @@ Use the transformed integrated API, not sampling at the light line. Method of images gives, with a=2 k1 d: - F_vertical = 1 + 3*(sin(a)-a*cos(a))/a^3 - F_parallel = 1 - (3/2)*(sin(a)/a + cos(a)/a^2 - sin(a)/a^3) +$$ +\begin{aligned} +F_\perp&=1+\frac{3(\sin a-a\cos a)}{a^3}\\ +F_\parallel&=1-\frac32\left(\frac{\sin a}{a}+\frac{\cos a}{a^2}-\frac{\sin a}{a^3}\right) +\end{aligned} +$$ At small a, the numerically stable limits are F_vertical=2-a^2/10+a^4/280+... and @@ -131,8 +145,9 @@ results detect wrong signs, missing prefactors and phase/unit mistakes. For a lossy half-space, the quasistatic leading evanescent corrections are: - Delta_vertical ~ 3/a^3 * Im[(epsilon2-epsilon1)/(epsilon2+epsilon1)] - Delta_parallel ~ (3/2)/a^3 * Im[(epsilon2-epsilon1)/(epsilon2+epsilon1)] +$$ +\Delta F_\perp\sim\frac{3}{a^3}\operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1},\qquad \Delta F_\parallel\sim\frac{3}{2a^3}\operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1} +$$ This is a mathematical asymptote, not a material/device-validity claim. diff --git a/docs/ENSEMBLE_RESOURCE_ESTIMATES.md b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md index 91430f5..5dc4b15 100644 --- a/docs/ENSEMBLE_RESOURCE_ESTIMATES.md +++ b/docs/ENSEMBLE_RESOURCE_ESTIMATES.md @@ -54,8 +54,12 @@ break that simplification. No finite-pixel EQE is claimed here. Let L be wavelength count, Z sampled depths, C thickness configurations, J finite layers, Q quadrature evaluations per orientation including refinements. - basis integrals = 2 * L * Z * C - layer-response evaluations = 2 * L * Z * C * Q * J +$$ +\begin{aligned} +N_{\operatorname{basis}}&=2LZC\\ +N_{\operatorname{layer}}&=2LZCQJ +\end{aligned} +$$ The estimator reports an adjustable arithmetic count, a streamed workspace, retained rate-array storage, and optional N-dipole direction storage/averaging @@ -136,10 +140,11 @@ factors. No multi-GPU capability is implied. `Hardware` optionally accepts USER-SUPPLIED sustained CPU/GPU GFLOP/s and GB/s. The estimator then reports only the idealized lower bound: - seconds >= max(assumed_FLOPs / sustained_FLOPs_per_second, - assumed_memory_traffic / sustained_bytes_per_second) +$$ +T\geq\max\left(\frac{\operatorname{FLOPs}_{\operatorname{assumed}}}{R_{\operatorname{FLOP}}},\frac{\operatorname{bytes}_{\operatorname{assumed}}}{R_{\operatorname{byte}}}\right) +$$ -This roofline lower bound is not elapsed-time prediction. If a grid does not +Here $R_{\operatorname{FLOP}}$ and $R_{\operatorname{byte}}$ are the user-supplied sustained rates. This roofline lower bound is not elapsed-time prediction. If a grid does not fit, the single-device bound does not mean the run is feasible. Without supplied rates, no CPU/GPU runtime is invented. The planning tool itself does not execute FDTD on either processor. diff --git a/docs/MATHEMATICAL_DERIVATIONS.md b/docs/MATHEMATICAL_DERIVATIONS.md index 42d3bfe..49bbde8 100644 --- a/docs/MATHEMATICAL_DERIVATIONS.md +++ b/docs/MATHEMATICAL_DERIVATIONS.md @@ -32,17 +32,13 @@ Assume linear, isotropic, nonmagnetic media and harmonic peak amplitudes with time dependence $$ -\mathbf E(\mathbf r,t)=\operatorname{Re} - \left[\mathbf E(\mathbf r)e^{-i\omega t}\right], -\qquad -k_0=\frac{2\pi}{\lambda_0}=\frac{\omega}{c}. +\mathbf E(\mathbf r,t)=\operatorname{Re} \left[\mathbf E(\mathbf r)e^{-i\omega t}\right], \qquad k_0=\frac{2\pi}{\lambda_0}=\frac{\omega}{c}. $$ The relative index and permittivity are $$ -N_j(\lambda_0)=n_j(\lambda_0)+i\kappa_j(\lambda_0), -\qquad \epsilon_j=N_j^2. +N_j(\lambda_0)=n_j(\lambda_0)+i\kappa_j(\lambda_0), \qquad \epsilon_j=N_j^2. $$ The symbol $\kappa$ denotes extinction coefficient, not a wavevector. @@ -57,8 +53,7 @@ formulas with $c$, $\epsilon_0$ and a dipole moment in C m. For a prescribed electric dipole $\mathbf p$, the homogeneous-host power is $$ -P_h=\frac{n_h\omega^4|\mathbf p|^2}{12\pi\epsilon_0c^3}. -\tag{1} +P_h=\frac{n_h\omega^4|\mathbf p|^2}{12\pi\epsilon_0c^3}. \tag{1} $$ This is a macroscopic lossless-host normalization for a **peak** amplitude, @@ -68,16 +63,10 @@ combining rates from different orientations or collection channels. The field radiated back onto the dipole changes source work. With the outgoing electric dyadic Green tensor defined by -$\mathbf E_{\rm ref}=\mu_0\omega^2\mathbf G_{\rm ref}\mathbf p$, +$\mathbf E_{\operatorname{ref}}=\mu_0\omega^2\mathbf G_{\operatorname{ref}}\mathbf p$, $$ -P_{\rm ref}=\frac{\omega}{2} - \operatorname{Im}\!\left[\mathbf p^*\cdot\mathbf E_{\rm ref}\right], -\qquad -F_{\rm EM}=1+\frac{6\pi}{k_h} - \operatorname{Im}\!\left[ - \hat{\mathbf e}^{\,*}\cdot\mathbf G_{\rm ref}\hat{\mathbf e}\right], -\tag{2} +P_{\operatorname{ref}}=\frac{\omega}{2} \operatorname{Im}\!\left[\mathbf p^*\cdot\mathbf E_{\operatorname{ref}}\right], \qquad F_{\operatorname{EM}}=1+\frac{6\pi}{k_h} \operatorname{Im}\!\left[ \hat{\mathbf e}^{\,*}\cdot\mathbf G_{\operatorname{ref}}\hat{\mathbf e}\right], \tag{2} $$ where $\hat{\mathbf e}=\mathbf p/|\mathbf p|$. The coefficient follows by @@ -96,12 +85,7 @@ Its squared magnitude is $1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2$. Normalizing its solid-angle integral gives $$ -\frac{1}{P_h}\frac{dP}{d\Omega} - =\frac{3}{8\pi} - \left(1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2\right), -\qquad -d\Omega=\sin\theta\,d\theta\,d\phi. -\tag{3} +\frac{1}{P_h}\frac{dP}{d\Omega} =\frac{3}{8\pi} \left(1-|\hat{\mathbf e}\cdot\hat{\mathbf r}|^2\right), \qquad d\Omega=\sin\theta\,d\theta\,d\phi. \tag{3} $$ Indeed, $\int d\Omega=4\pi$ and @@ -113,12 +97,7 @@ Project the dipole onto the orthogonal far-field polarization axes before squaring: $$ -\frac{dP_s}{P_h\,d\Omega}=\frac{3}{8\pi} - |\hat{\mathbf e}\cdot\hat{\boldsymbol\phi}|^2, -\qquad -\frac{dP_p}{P_h\,d\Omega}=\frac{3}{8\pi} - |\hat{\mathbf e}\cdot\hat{\boldsymbol\theta}|^2. -\tag{4} +\frac{dP_s}{P_h\,d\Omega}=\frac{3}{8\pi} |\hat{\mathbf e}\cdot\hat{\boldsymbol\phi}|^2, \qquad \frac{dP_p}{P_h\,d\Omega}=\frac{3}{8\pi} |\hat{\mathbf e}\cdot\hat{\boldsymbol\theta}|^2. \tag{4} $$ | Dipole distribution | $s$ factor, excluding $3/(8\pi)$ | $p$ factor, excluding $3/(8\pi)$ | @@ -139,11 +118,7 @@ with $\epsilon_1=n_h^2$. Medium 2 occupies $z<0$. For each in-plane wavevector magnitude, define $$ -q=\frac{k_\parallel}{k_h},\qquad -w=\sqrt{1-q^2},\qquad -v=\sqrt{\frac{\epsilon_2}{\epsilon_1}-q^2}, -\qquad a=2k_hd=\frac{4\pi n_hd}{\lambda_0}. -\tag{5} +q=\frac{k_\parallel}{k_h},\qquad w=\sqrt{1-q^2},\qquad v=\sqrt{\frac{\epsilon_2}{\epsilon_1}-q^2}, \qquad a=2k_hd=\frac{4\pi n_hd}{\lambda_0}. \tag{5} $$ On the real $q$ axis, $w$ is positive real below the host light line @@ -152,10 +127,7 @@ outgoing/decaying prescription. Electric-field Fresnel reflection amplitudes are $$ -r_s=\frac{w-v}{w+v},\qquad -r_p=\frac{\epsilon_2w-\epsilon_1v} - {\epsilon_2w+\epsilon_1v}. -\tag{6} +r_s=\frac{w-v}{w+v},\qquad r_p=\frac{\epsilon_2w-\epsilon_1v} {\epsilon_2w+\epsilon_1v}. \tag{6} $$ Identical media give $r_s=r_p=0$. The explicit perfect-electric-conductor @@ -166,15 +138,11 @@ The angular-spectrum representation of the reflected Green tensor, followed by azimuthal integration and division by (1), gives $$ -F_\perp=1+\frac32\operatorname{Re}\int_0^\infty - \frac{q^3}{w}r_p e^{iaw}\,dq, -\tag{7} +F_\perp=1+\frac32\operatorname{Re}\int_0^\infty \frac{q^3}{w}r_p e^{iaw}\,dq, \tag{7} $$ $$ -F_\parallel=1+\frac34\operatorname{Re}\int_0^\infty - \frac{q}{w}\left(r_s-w^2r_p\right)e^{iaw}\,dq. -\tag{8} +F_\parallel=1+\frac34\operatorname{Re}\int_0^\infty \frac{q}{w}\left(r_s-w^2r_p\right)e^{iaw}\,dq. \tag{8} $$ Here $\perp$ means dipole normal to the interface, not the $s$ polarization. @@ -186,7 +154,7 @@ The ideal-mirror tests below independently constrain this sign. The $q<1$ / $q>1$ split classifies propagation **in the host**. A host-evanescent component can propagate in a higher-index substrate. -Thus neither the evanescent correction nor $F_{\rm EM}$ alone is an +Thus neither the evanescent correction nor $F_{\operatorname{EM}}$ alone is an absorption fraction or an external collection efficiency. ### 4.1 Removing the host light-line singularity @@ -196,38 +164,28 @@ sampling. In the propagating sector use $u=\sqrt{1-q^2}$, so $q\,dq=-u\,du$. Reversing the integration limits yields $$ -\Delta F_{\perp,\rm prop}=\frac32\int_0^1 - (1-u^2)\operatorname{Re}[r_pe^{iau}]\,du, -\tag{9} +\Delta F_{\perp,\operatorname{prop}}=\frac32\int_0^1 (1-u^2)\operatorname{Re}[r_pe^{iau}]\,du, \tag{9} $$ $$ -\Delta F_{\parallel,\rm prop}=\frac34\int_0^1 - \operatorname{Re}[(r_s-u^2r_p)e^{iau}]\,du. -\tag{10} +\Delta F_{\parallel,\operatorname{prop}}=\frac34\int_0^1 \operatorname{Re}[(r_s-u^2r_p)e^{iau}]\,du. \tag{10} $$ For the evanescent sector set $t=\sqrt{q^2-1}$, then $w=it$, $q\,dq=t\,dt$ and $e^{iaw}=e^{-at}$. Finally let $x=at$: $$ -\Delta F_{\perp,\rm evan}=\frac{3}{2a}\int_0^\infty - \left[1+\left(\frac{x}{a}\right)^2\right] - e^{-x}\operatorname{Im}r_p\,dx, -\tag{11} +\Delta F_{\perp,\operatorname{evan}}=\frac{3}{2a}\int_0^\infty \left[1+\left(\frac{x}{a}\right)^2\right] e^{-x}\operatorname{Im}r_p\,dx, \tag{11} $$ $$ -\Delta F_{\parallel,\rm evan}=\frac{3}{4a}\int_0^\infty - e^{-x}\operatorname{Im}\left[ - r_s+\left(\frac{x}{a}\right)^2r_p\right]dx. -\tag{12} +\Delta F_{\parallel,\operatorname{evan}}=\frac{3}{4a}\int_0^\infty e^{-x}\operatorname{Im}\left[ r_s+\left(\frac{x}{a}\right)^2r_p\right]dx. \tag{12} $$ The coefficients in each transformed integral are evaluated at that integral's $q$. The Jacobians cancel the host light-line singularity; scaling by $a$ also makes the dominant decay envelope $e^{-x}$ independent of separation. -The full answer is $1+\Delta F_{\rm prop}+\Delta F_{\rm evan}$. +The full answer is $1+\Delta F_{\operatorname{prop}}+\Delta F_{\operatorname{evan}}$. ### 4.2 Branches, resonances and error gates @@ -235,9 +193,7 @@ Write $\rho=\epsilon_2/\epsilon_1$. For a lossless dielectric, its light line appears at $$ -u_b=\sqrt{1-\rho}\quad(0<\rho<1), -\qquad x_b=a\sqrt{\rho-1}\quad(\rho>1). -\tag{13} +u_b=\sqrt{1-\rho}\quad(0<\rho<1), \qquad x_b=a\sqrt{\rho-1}\quad(\rho>1). \tag{13} $$ Splitting there helps resolve changes in the substrate root. Weak loss moves @@ -245,16 +201,15 @@ these branch locations off the real axis; the integration intervals are split around their projected centers and widths. A $p$ surface-mode denominator in (6) gives -$q_{\rm sp}^2=\rho/(\rho+1)$ and hence +$q_{\operatorname{sp}}^2=\rho/(\rho+1)$ and hence $$ -t_{\rm sp}=\sqrt{-\frac{1}{\rho+1}},\qquad x_{\rm sp}=at_{\rm sp}. -\tag{14} +t_{\operatorname{sp}}=\sqrt{-\frac{1}{\rho+1}},\qquad x_{\operatorname{sp}}=at_{\operatorname{sp}}. \tag{14} $$ For supported lossy-metal cases, intervals around -$\operatorname{Re}x_{\rm sp}$ span geometrically increasing multiples of -$|\operatorname{Im}x_{\rm sp}|$. This resolves narrow resonances and wings +$\operatorname{Re}x_{\operatorname{sp}}$ span geometrically increasing multiples of +$|\operatorname{Im}x_{\operatorname{sp}}|$. This resolves narrow resonances and wings better than a single large adaptive interval. The implementation integrates to infinity, using an initial tail split @@ -264,11 +219,7 @@ Per-interval absolute budgets share the scaled absolute tolerance. For each orientation it checks the summed quadrature error estimate and $$ -\delta_{\rm refine}= - |\Delta F_{\rm prop}^{(2)}-\Delta F_{\rm prop}^{(1)}| - +|\Delta F_{\rm evan}^{(2)}-\Delta F_{\rm evan}^{(1)}| - \leq \max(\mathrm{atol},\mathrm{rtol}|F|). -\tag{15} +\delta_{\operatorname{refine}}= |\Delta F_{\operatorname{prop}}^{(2)}-\Delta F_{\operatorname{prop}}^{(1)}| +|\Delta F_{\operatorname{evan}}^{(2)}-\Delta F_{\operatorname{evan}}^{(1)}| \leq \max(\operatorname{atol},\operatorname{rtol}|F|). \tag{15} $$ Failures, nonfinite results and negative total rates raise errors; rates are @@ -287,20 +238,13 @@ that local continuum material constants remain valid at molecular distances. For an ideal mirror, image-dipole interference gives $$ -F_\perp=1+3\frac{\sin a-a\cos a}{a^3}, -\qquad -F_\parallel=1-\frac32\left( - \frac{\sin a}{a}+\frac{\cos a}{a^2}-\frac{\sin a}{a^3}\right). -\tag{16} +F_\perp=1+3\frac{\sin a-a\cos a}{a^3}, \qquad F_\parallel=1-\frac32\left( \frac{\sin a}{a}+\frac{\cos a}{a^2}-\frac{\sin a}{a^3}\right). \tag{16} $$ Series expansion avoids cancellation near the surface: $$ -F_\perp=2-\frac{a^2}{10}+\frac{a^4}{280}+O(a^6), -\qquad -F_\parallel=\frac{a^2}{5}-\frac{3a^4}{280}+O(a^6). -\tag{17} +F_\perp=2-\frac{a^2}{10}+\frac{a^4}{280}+O(a^6), \qquad F_\parallel=\frac{a^2}{5}-\frac{3a^4}{280}+O(a^6). \tag{17} $$ The vertical image reinforces the source while the horizontal image cancels @@ -312,12 +256,7 @@ $r_s\to0$, $r_p\to(\epsilon_2-\epsilon_1)/(\epsilon_2+\epsilon_1)$. Using $\int_0^\infty x^2e^{-x}dx=2$ in (11)-(12) yields $$ -\Delta F_\perp\sim\frac{3}{a^3} - \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}, -\qquad -\Delta F_\parallel\sim\frac{3}{2a^3} - \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}. -\tag{18} +\Delta F_\perp\sim\frac{3}{a^3} \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}, \qquad \Delta F_\parallel\sim\frac{3}{2a^3} \operatorname{Im}\frac{\epsilon_2-\epsilon_1}{\epsilon_2+\epsilon_1}. \tag{18} $$ This explains the leading inverse-cube dependence and two-to-one orientation @@ -330,14 +269,11 @@ For tabulated data between $\lambda_i$ and $\lambda_{i+1}$, the current interpolator applies independently to $n$ and $\kappa$: $$ -f(\lambda)=f_i+ - \frac{\lambda-\lambda_i}{\lambda_{i+1}-\lambda_i}(f_{i+1}-f_i), -\qquad f\in\{n,\kappa\}. -\tag{19} +f(\lambda)=f_i+ \frac{\lambda-\lambda_i}{\lambda_{i+1}-\lambda_i}(f_{i+1}-f_i), \qquad f\in\{n,\kappa\}. \tag{19} $$ Native database wavelengths in micrometers are converted to nm: -$\lambda_{\rm nm}=1000\lambda_{\rm um}$. +$\lambda_{\operatorname{nm}}=1000\lambda_{\operatorname{um}}$. No extrapolation is allowed. Separate $n$ and $\kappa$ tables must both cover the query, so their usable domain is their intersection. Missing extinction data require an explicit, recorded lossless assumption; missing data are not @@ -347,12 +283,8 @@ The supported native Sellmeier records use $\lambda$ in **micrometers**: $$ \begin{aligned} -\text{formula 1:}\quad -n^2(\lambda)=1+c_0+\sum_m - \frac{B_m\lambda^2}{\lambda^2-C_m^2},\\ -\text{formula 2:}\quad -n^2(\lambda)=1+c_0+\sum_m - \frac{B_m\lambda^2}{\lambda^2-C_m}. +\text{formula 1:}\; n^2(\lambda)=1+c_0+\sum_m \frac{B_m\lambda^2}{\lambda^2-C_m^2}, \\ +\text{formula 2:}\; n^2(\lambda)=1+c_0+\sum_m \frac{B_m\lambda^2}{\lambda^2-C_m}. \end{aligned} \tag{20} $$ @@ -365,9 +297,7 @@ Each layer has its own $N_j(\lambda)$ and thickness $d_j$. For a plane-wave component with conserved $k_\parallel$, $$ -k_{z,j}=\sqrt{N_j^2k_0^2-k_\parallel^2}, -\qquad \delta_j=k_{z,j}d_j. -\tag{21} +k_{z,j}=\sqrt{N_j^2k_0^2-k_\parallel^2}, \qquad \delta_j=k_{z,j}d_j. \tag{21} $$ Thickness changes phase and attenuation, while dispersion changes propagation @@ -378,9 +308,7 @@ The implemented diagnostic calls the coherent transfer-matrix calculation for both incident polarizations and averages their **powers**: $$ -R=\frac{R_s+R_p}{2},\qquad -T=\frac{T_s+T_p}{2},\qquad A=1-R-T. -\tag{22} +R=\frac{R_s+R_p}{2},\qquad T=\frac{T_s+T_p}{2},\qquad A=1-R-T. \tag{22} $$ For the enforced lossless exterior media, power transmission includes @@ -404,34 +332,23 @@ Define an intrinsic radiative rate in the homogeneous reference, $\Gamma_{r,h}$, and a fixed intrinsic nonradiative rate $\Gamma_{nr}$: $$ -q_0=\frac{\Gamma_{r,h}}{\Gamma_{r,h}+\Gamma_{nr}}, -\qquad -\Gamma_{nr}=\Gamma_{r,h}\frac{1-q_0}{q_0} -\quad(0