update

examples/calc-size-averaged-Qext.py

@@ -33,7 +33,7 @@ from optical_constants import read_refractive_index_from_yaml as get_index
 def gauss(x, mu, sigma):
     return 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (x - mu)**2 / (2 * sigma**2))
 
-fill_factor = 0.8
+# fill_factor = 0.8
 
 from_WL = 300
 to_WL = 1100
@@ -74,36 +74,36 @@ index_TiO2 = get_index("refractiveindex_info/TiO2-Sarkar.yml", WLs, units='nm')
 # index_TiO2[:,1] += 0.5j
 # index_Au[:,1] = index_Au[:,1]* fill_factor + (1-fill_factor)
 
-print("Au index before correction, max = ", np.max(np.imag(index_Au[:,1])))
-# Taking into account gold free electrons damping
-# contributed by surface scattering and bulk dumping
-# See eq 1 in [1] -> doi: https://doi.org/10.1186/s11671-018-2670-7
-eps_exp = index_Au[:,1]**2
-c=299792458  # speed of light
-h= 4.135667516e-15  # eV*s, Planck constant
-w = h*c/(WLs*1e-9)  # eV, frequency
-w_p =8.55  # eV, gold plasmon frequency
-g_b = 18.4e-3  # eV, bulk dumping
-v_f = 1.4e6  # m/s, Fermi velocity of electrons in gold
-A = 1.33  # fit parameter for 16 nm gold shell, see Table 2 in [1]
-L_b = (4.*((core_r+inner_shell_h)**3 - core_r**3)/
-       (3.*((core_r+inner_shell_h)**2 + core_r**2)))
-# g_s = v_f/L_b  # eq 2 in [1]
-g_s = h*A*v_f/(inner_shell_h*1e-9)  # eq 4 in [1]
-# g_b *= 0.2
-# g_s *= 0.5
-eps_Au = (eps_exp
-          +
-          w_p**2 / ( w * (w + 1j*g_b) )
-          -
-          w_p**2 / ( w * (w + 1j*(g_b+g_s)) )
-          )
+# print("Au index before correction, max = ", np.max(np.imag(index_Au[:,1])))
+# # Taking into account gold free electrons damping
+# # contributed by surface scattering and bulk dumping
+# # See eq 1 in [1] -> doi: https://doi.org/10.1186/s11671-018-2670-7
+# eps_exp = index_Au[:,1]**2
+# c=299792458  # speed of light
+# h= 4.135667516e-15  # eV*s, Planck constant
+# w = h*c/(WLs*1e-9)  # eV, frequency
+# w_p =8.55  # eV, gold plasmon frequency
+# g_b = 18.4e-3  # eV, bulk dumping
+# v_f = 1.4e6  # m/s, Fermi velocity of electrons in gold
+# A = 1.33  # fit parameter for 16 nm gold shell, see Table 2 in [1]
+# L_b = (4.*((core_r+inner_shell_h)**3 - core_r**3)/
+#        (3.*((core_r+inner_shell_h)**2 + core_r**2)))
+# # g_s = v_f/L_b  # eq 2 in [1]
+# g_s = h*A*v_f/(inner_shell_h*1e-9)  # eq 4 in [1]
+# # g_b *= 0.2
+# # g_s *= 0.5
+# eps_Au = (eps_exp
+#           +
+#           w_p**2 / ( w * (w + 1j*g_b) )
+#           -
+#           w_p**2 / ( w * (w + 1j*(g_b+g_s)) )
+#           )
 # index_Au[:,1] = np.sqrt(eps_Au)
 index_Au[:,1] += 1.6j
-print(f"L_b={L_b}")
-# print(w)
-print(f"g_s={g_s}, g_b={g_b} ")
-print("Au index after, max = ", np.max(np.imag(index_Au[:,1])))
+# print(f"L_b={L_b}")
+# # print(w)
+# print(f"g_s={g_s}, g_b={g_b} ")
+# print("Au index after, max = ", np.max(np.imag(index_Au[:,1])))
 x = np.ones((3), dtype = np.float64)
 m = np.ones((3), dtype = np.complex128)
 

examples/calc-spectra.py

@@ -46,7 +46,7 @@ x = np.ones((1), dtype = np.float64)
 m = np.ones((1), dtype = np.complex128)
 
 
-core_r = 45000
+core_r = 100
 
 
 Qsca_vec = []

examples/optical_constants.py

@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+# File from Smuthi project https://gitlab.com/AmosEgel/smuthi/-/blob/master/smuthi/utility/optical_constants.py
+"""Provide functionality to read optical constants in format provided by `refractiveindex.info <https://refractiveindex.info/>`_ website"""
+
+from scipy.interpolate import interp1d
+import numpy as np
+import yaml
+
+def read_refractive_index_from_yaml(filename, vacuum_wavelength, units="mkm", kind=1):
+    """Read optical constants in format provided by refractiveindex.info website.
+
+    Args:
+            filename (str): path and file name for yaml data
+                            downloaded from refractiveindex.info
+            vacuum_wavelength (float or np.array): wavelengths where refractive
+                                                   index data is needed
+            units (str): units for wavelength. currently, microns ('mkm' or 'um')
+                         and nanometers ('nm') can be selected
+            kind (int): order of interpolation
+
+    Returns:
+            A pair (or np.array of pairs) of wavelength and
+            corresponding refractive index (complex)
+    """
+    if units == "nm":
+        factor = 1000
+    elif units in ("mkm", "um"):
+        factor = 1
+    else:
+        raise NotImplementedError("Converting wavelength into '"+units
+                                  +"' units for refractive index data"
+                                  +" was not implemented.")
+
+    with open(filename) as f:
+        the_file = yaml.load(f)['DATA'][0]
+    data_type = the_file['type']
+    if data_type != 'tabulated nk':
+        raise NotImplementedError("Input data type '"+data_type
+                                  +"' available in file "+filename
+                                  +" was not implemented.")
+    data = the_file['data'].splitlines()
+    data_split = []
+    for wl in data:
+        data_split.append(wl.split())
+    data_num = []
+    for wl in data_split:
+        record = []
+        for val in wl:
+            record.append(float(val))
+        data_num.append(record)
+    data_np = np.array(data_num)
+    data_wl = data_np[:,0]*factor
+    index_re = data_np[:,1]
+    index_im = data_np[:,2]
+    f_re = interp1d(data_wl, index_re, kind=kind)
+    f_im = interp1d(data_wl, index_im, kind=kind)
+    data_out = np.transpose(np.vstack((vacuum_wavelength, f_re(vacuum_wavelength)+f_im(vacuum_wavelength)*1j)))
+    if len(data_out) == 1:
+        return data_out[0]
+    return data_out

tests/mpmath_riccati_bessel.py

@@ -16,7 +16,7 @@ def psi(n,z):
 # Riccati-Bessel -z*y_n(z)
 def xi(n,z):
     return -mp.sqrt( (mp.pi * z)/2 ) * mp.autoprec(mp.bessely)(n+1/2,z)
-# Riccati-Bessel psi - i* xi
+# Riccati-Bessel psi - i* xi        z*j_n(z) + i* z*j_n(z)
 def ksi(n,z):
     return psi(n,z) - 1.j * xi(n,z)
 

vue-cli3-webapp/src/App.vue

@@ -124,6 +124,7 @@
   // ).then(bytes =>
   //         console.log(bytes)
   // );
+  import { mapState } from 'vuex'
 
   import nmiejs from './nmiejs.js';
   // Test nmiejs if working
@@ -180,6 +181,10 @@
       ReactiveChart,
       ShowInfo
     },
+    computed:
+        mapState([
+          'simulationSetup'
+        ]),
     data() {
         return {
           nmie: undefined,
@@ -194,23 +199,6 @@
           source_units_prev: 'nm',
           materials: [],
           isSourceOtherUnits: false,
-          simulationSetup: {
-            hostIndex: 1,
-            stepWL: 2,
-            fromWL: "100*(1+2)",
-            toWL: 1000,
-            layers: [
-              {
-                R: 100.0,
-                material: 'nk',
-                isMaterialLoaded:true,
-                reN: 4.0,
-                imN: 0.01,
-                index: 0
-              }
-            ],
-            total_mode_n: 4
-          },
           simulationRuntime: {
             r_units: 'nm',
             r_source_units: 'nm',
@@ -335,6 +323,28 @@
 
     },
     watch: {
+      'simulationSetup.layers': {
+        handler: function () {
+          this.checkConflictsWL();
+        },
+        deep: true
+      },
+      'simulationSetup.toWL': {
+        handler: function () {
+          for (let i = 0; i < this.simulationSetup.layers.length; i++) {
+            this.simulationSetup.layers[i].isMaterial_hasConflict = false;
+          }
+          this.checkConflictsWL();
+        }
+      },
+      'simulationSetup.fromWL': {
+        handler: function () {
+          for (let i = 0; i < this.simulationSetup.layers.length; i++) {
+            this.simulationSetup.layers[i].isMaterial_hasConflict = false;
+          }
+          this.checkConflictsWL();
+        }
+      },
       plotSelector: {
           handler: function () {
             this.plotResults();
@@ -420,6 +430,29 @@
       }
     },
     methods: {
+      checkConflictsWL (){
+        this.simulationSetup.fromWL_hasConflict = false;
+        this.simulationSetup.toWL_hasConflict = false;
+        let fromWL = this.convertUnits2nm(this.source_units, this.simulationSetup.fromWL);
+        let toWL = this.convertUnits2nm(this.source_units, this.simulationSetup.toWL);
+
+        for (let i = 0; i < this.simulationSetup.layers.length; i++) {
+          // console.log('checking WL conflict layer', i);
+          if (this.simulationSetup.layers[i].material=='nk') continue;
+          // console.log(fromWL, ' and ', this.simulationSetup.layers[i].spline_n.xs[0]);
+          // console.log(toWL, ' and ', this.simulationSetup.layers[i]
+          //     .spline_n.xs[this.simulationSetup.layers[i].spline_n.xs.length-1] );
+          if (fromWL < this.simulationSetup.layers[i].spline_n.xs[0]) {
+            this.simulationSetup.fromWL_hasConflict = true;
+            this.simulationSetup.layers[i].isMaterial_hasConflict = true;
+          }
+          if ( toWL > this.simulationSetup.layers[i]
+              .spline_n.xs[this.simulationSetup.layers[i].spline_n.xs.length-1] ) {
+            this.simulationSetup.toWL_hasConflict = true;
+            this.simulationSetup.layers[i].isMaterial_hasConflict = true;
+          }
+        }
+      },
       handleResize() {
           this.window.width = window.innerWidth;
           this.window.height = window.innerHeight*0.8;
@@ -485,7 +518,34 @@
         });
       },
       runSimulation: function() {
-          this.$buefy.notification.open({
+          if (this.simulationSetup.toWL_hasConflict ||
+          this.simulationSetup.fromWL_hasConflict) {
+            this.$buefy.notification.open({
+              duration: 2000,
+              message: 'Simulation has input conflicts, fixing...',
+              type: 'is-danger',
+              position: 'is-bottom-left',
+            });
+          }
+
+        for (let i = 0; i < this.simulationSetup.layers.length; i++) {
+          if (this.simulationSetup.layers[i].material=='nk') continue;
+          let fromWLnm = this.convertUnits2nm(this.source_units, this.simulationSetup.fromWL);
+          let toWLnm = this.convertUnits2nm(this.source_units, this.simulationSetup.toWL);
+          let fromWL_material = this.simulationSetup.layers[i].spline_n.xs[0];
+          if (fromWLnm < fromWL_material ) {
+            this.simulationSetup.fromWL = this.convertUnitsFrom_nm(
+                this.source_units, fromWL_material);
+          }
+          let toWL_material = this.simulationSetup.layers[i]
+              .spline_n.xs[this.simulationSetup.layers[i].spline_n.xs.length-1];
+          if ( toWLnm > toWL_material  ) {
+            this.simulationSetup.toWL = this.convertUnitsFrom_nm(
+                this.source_units, toWL_material);
+          }
+        }
+
+        this.$buefy.notification.open({
             duration: 200,
             message: 'Simulation was started!',
             type: 'is-danger',
@@ -567,6 +627,9 @@
               if (layer.spline_n.xs[0] > WL
                       || layer.spline_n.xs[layer.spline_n.xs.length-1] < WL) {
                 this.notifyDanger('ERROR!!! Source parameters are out of material spectral range!');
+                layer.isMaterial_hasConflict = true;
+                this.simulationSetup.toWL_hasConflict = true;
+                this.simulationSetup.fromWL_hasConflict = true;
                 return;
               }
               layer.reN = layer.spline_n.at(WL);

vue-cli3-webapp/src/components/GetLayerParameters.vue

@@ -1,5 +1,5 @@
 <template>
-    <div class="field">
+    <div class="field" >
         <div class="field is-horizontal layer">
             <div class="field-label is-normal">
                 <label class="label lnorm">
@@ -28,15 +28,19 @@
                         <span class="tooltiptext">Layer thickness</span>
 
                     </div>
-                    <b-select v-model="layer.material">
-                        <option v-if="index == 0" value="PEC" disabled> PEC</option>
-                        <option value="nk">nk-constant</option>
-                        <option v-for="material in filteredMaterials" v-bind:key="material.name"
-                                v-bind:value="material.name">
-                            {{material.name}} ({{material.spline_n.xs[0]}} -
-                                               {{material.spline_n.xs[material.spline_n.xs.length-1]}} nm)
-                        </option>
-                    </b-select>
+                    <b-field :type="{ 'is-danger': layer.isMaterial_hasConflict }"
+                             :message="{ 'Input conflict': layer.isMaterial_hasConflict }"
+                             style="margin-bottom: 0.5rem">
+                      <b-select v-model="layer.material">
+                            <option v-if="index == 0" value="PEC" disabled> PEC</option>
+                            <option value="nk">nk-constant</option>
+                            <option v-for="material in filteredMaterials" v-bind:key="material.name"
+                                    v-bind:value="material.name">
+                                {{material.name}} ({{material.spline_n.xs[0]}} -
+                                                   {{material.spline_n.xs[material.spline_n.xs.length-1]}} nm)
+                            </option>
+                      </b-select>
+                    </b-field>
                     <div v-if="layer.isMaterialLoaded">
                         <font-awesome-icon icon="check-circle" class="has-text-success"/>
                     </div>
@@ -91,28 +95,29 @@
             }
         },
         watch: {
-            'layer.material': {
-                handler: function () {
-                    // console.log('update');
-                    this.layer.isMaterialLoaded = false;
-                    if (this.layer.material == 'nk') {
-                        this.isDisabled = false;
-                        this.layer.reN = 4;
-                        this.layer.imN = 0.01;
-                        this.layer.isMaterialLoaded = true;
-                    } else {
-                        this.isDisabled = true;
-                    }
-                    for (const mat of this.materials) {
-                        if (mat.name != this.layer.material) continue;
-                        this.layer.spline_n = mat.spline_n;
-                        this.layer.spline_k = mat.spline_k;
-                        console.log(mat.spline_n.at(500));
-                        console.log(this.layer.spline_n.at(500));
-                        this.layer.isMaterialLoaded = true;
-                    }
-                }
+          'layer.material': {
+            handler: function () {
+              // console.log('update');
+              this.layer.isMaterialLoaded = false;
+              if (this.layer.material == 'nk') {
+                this.isDisabled = false;
+                this.layer.reN = 4;
+                this.layer.imN = 0.01;
+                this.layer.isMaterialLoaded = true;
+              } else {
+                this.isDisabled = true;
+              }
+              for (const mat of this.materials) {
+                if (mat.name != this.layer.material) continue;
+                this.layer.spline_n = mat.spline_n;
+                this.layer.spline_k = mat.spline_k;
+                console.log(mat.spline_n.at(500));
+                console.log(this.layer.spline_n.at(500));
+                this.layer.isMaterialLoaded = true;
+              }
+              this.layer.isMaterial_hasConflict = false;
             }
+          }
         },
         props: ['layer', 'particle', 'index', 'units', 'materials']
     }
@@ -148,7 +153,7 @@
 
         /* Position the tooltip text - see examples below! */
         position: absolute;
-        z-index: 1;
+        z-index: 4;
 
         bottom: 120%;
         left: 0%;
@@ -159,6 +164,6 @@
         visibility: visible;
     }
     .nk-input {
-        margin-bottom: 12px;
+        margin-bottom: 0.25rem;
     }
 </style>

vue-cli3-webapp/src/components/GetMaterials.vue

@@ -179,6 +179,12 @@
                 if (mat.name.includes('Johnson')) mat.isUsed=false;
                 if (mat.name.includes('Aspnes')) mat.isUsed=false;
             }
+          try { // wake-up proxcors
+            let tryURL = 'https://proxcors.herokuapp.com/https://refractiveindex.info/database/data/main/Ag/McPeak.yml';
+            fetch(tryURL);
+          } catch (e) {
+            console.log(e);
+          }
         },
         watch: {
             materials: {
@@ -329,7 +335,7 @@
                 let Ag_data;
 
                 try {
-                    let proxyUrl = 'https://cors-anywhere.herokuapp.com/';
+                    let proxyUrl = 'https://proxcors.herokuapp.com/';
                     if (URL.includes('https://refractiveindex.info')) URL = proxyUrl+URL;
                     let response = await fetch(URL);
                     let Ag_data = await response.text();

vue-cli3-webapp/src/components/GetParticleParameters.vue

@@ -32,6 +32,7 @@
 
 <script>
     import GetLayerParameters from "./GetLayerParameters.vue";
+    import { mapState } from 'vuex'
     export default {
         name: "GetParticleParameters",
         components: {
@@ -45,7 +46,11 @@
                 index: 1
             }
         },
-        watch: {
+      computed:
+          mapState([
+            'simulationSetup'
+          ]),
+      watch: {
             // emit updated values
             particle: {
                 handler: function () {
@@ -64,18 +69,20 @@
             // },
             layersNum: {
                 handler: function () {
-                    while (this.layersNum < this.layers.length) {
-                        this.layers.pop();
+                    while (this.layersNum < this.simulationSetup.layers.length) {
+                        this.simulationSetup.layers.pop();
                     }
-                    let r_prev = this.layers[0].R;
-                    while (this.layersNum > this.layers.length) {
+                    let r_prev = this.simulationSetup.layers[0].R;
+                    while (this.layersNum > this.simulationSetup.layers.length) {
                         // r_prev = r_prev*1.1;
-                        this.layers.push(
+                        this.simulationSetup.layers.push(
                             {
                                 R: r_prev*0.1,
                                 material: 'nk',
                                 isMaterialLoaded:true,
-                                reN: 4.0,
+                              isMaterial_hasConflict:false,
+
+                              reN: 4.0,
                                 imN: 0.01,
                                 index: 1,
                                 spline_n: undefined,
@@ -83,8 +90,8 @@
                             }
                         );
                     }
-                    for (let i = 0; i < this.layers.length; i++) {
-                        this.layers[i].index = i;
+                    for (let i = 0; i < this.simulationSetup.layers.length; i++) {
+                        this.simulationSetup.layers[i].index = i;
                     }
                 }
             },

vue-cli3-webapp/src/components/GetSourceParameters.vue

@@ -12,9 +12,11 @@
             <div class="field is-grouped is-grouped-multiline">
                 <input-with-units title="from" v-bind:units="source_units"
                                   v-bind:value="fromWLLocal"
+                                  v-bind:has-conflict="simulationSetup.fromWL_hasConflict"
                                   @newdata="fromWLLocal=$event"/>
                 <input-with-units title="to" v-bind:units="source_units"
                                   v-bind:value="toWLLocal"
+                                  v-bind:has-conflict="simulationSetup.toWL_hasConflict"
                                   @newdata="toWLLocal=$event"/>
                 <div>
                 <div v-if="isStep">
@@ -36,6 +38,7 @@
 
 <script>
 import InputWithUnits from "./InputWithUnits.vue";
+import { mapState } from 'vuex'
 export default {
   name: "GetSourceParameters",
   components: {
@@ -54,6 +57,10 @@ export default {
   methods: {
 
   },
+  computed:
+    mapState([
+    'simulationSetup'
+  ]),
   watch: {
     isStep: {
       handler: function () {

vue-cli3-webapp/src/components/InputWithUnits.vue

@@ -1,5 +1,5 @@
 <template>
-    <div class="field has-addons">
+    <div class="field has-addons" style="margin-bottom: 0.25rem">
         <p class="control">
             <a class="button is-static input-with-units-title">
                 {{ title }}
@@ -14,12 +14,15 @@
         </div>
         <div v-else>
             <p class="control">
+              <b-field :type="{ 'is-danger': hasConflict }"
+                       :message="{ 'Input conflict': hasConflict }">
                 <b-autocomplete v-on:focus="handleFocus($event)" v-on:blur="handleBlur"
                          v-on:keyup.native.enter="focusNext($event)" v-model="showValue"
                                 :data="expr_list"
                                 :open-on-focus="true"
                                 class="input-with-units-value"/>
 
+              </b-field>
             </p>
         </div>
         <p class="control">
@@ -102,6 +105,11 @@
                     this.isDisabledLocal = this.isDisabled;
                 }
             },
+          hasConflict: {
+            handler: function () {
+              this.hasConflictLocal = this.hasConflict;
+            }
+          },
           value: {
             handler: function () {
               this.showValue = this.value;
@@ -117,14 +125,15 @@
               expr_list: [],
               showValue: math.evaluate(this.value),
               isFocused: false,
-                isDisabledLocal: false,
+              isDisabledLocal: false,
+              hasConflictLocal: this.hasConflict
             }
         },
       mounted() {
         this.$emit('newdata',math.evaluate(this.value));
 
       },
-        props: ['title', 'units', 'value', 'isDisabled']
+        props: ['title', 'units', 'value', 'isDisabled', 'hasConflict']
     }
 </script>
 
@@ -136,6 +145,7 @@
 
 .input-with-units-value {
     width:10rem;
+   z-index: 2;
 }
 .input-with-units-units {
     width:3rem;

vue-cli3-webapp/src/store.js

@@ -5,7 +5,26 @@ Vue.use(Vuex)
 
 export default new Vuex.Store({
     state: {
-        count: 0
+        simulationSetup: {
+            hostIndex: 1,
+            stepWL: 2,
+            fromWL: "100*(1+2)",
+            fromWL_hasConflict: false,
+            toWL: 1000,
+            toWL_hasConflict:false,
+            layers: [
+                {
+                    R: 100.0,
+                    material: 'nk',
+                    isMaterialLoaded:true,
+                    isMaterial_hasConflict:false,
+                    reN: "sqrt(16)",
+                    imN: 0.01,
+                    index: 0
+                }
+            ],
+            total_mode_n: 4
+        },
     },
     mutations: {
         increment (state) {