add RunFieldCalculationCartesian method to nmie.hpp and link it to JS module

CMakeLists.txt

@@ -21,11 +21,6 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 add_compile_options(-W -Wall -pedantic -Werror)
 add_compile_options(-funroll-loops -fstrict-aliasing)
 
-# Set options
-option(ENABLE_MP "Use multiple precision" OFF)
-if (${ENABLE_MP})
-    add_compile_options(-DMULTI_PRECISION=100)
-endif ()
 
 # compiler details
 message("  C++ Compiler: ${CMAKE_CXX_COMPILER_ID} "
@@ -42,17 +37,24 @@ if (USE_STATIC_LIBRARIES)
 endif ()
 set(Boost_USE_MULTITHREADED OFF)
 set(Boost_USE_STATIC_RUNTIME OFF)
-find_package(Boost REQUIRED)
-if (Boost_INCLUDE_DIRS)
-    if (${Boost_VERSION} VERSION_LESS 1.60.0)
-        message(FATAL_ERROR
-                "Found Boost library is too old; required is version 1.60.0 or newer!")
+find_package(Boost)
+# Set options
+option(ENABLE_MP "Use multiple precision" OFF)
+if (Boost_FOUND)
+    if (${ENABLE_MP})
+        add_compile_options(-DMULTI_PRECISION=100)
     endif ()
-    message("Found Boost include dir: ${Boost_INCLUDE_DIR}")
-    message("Found Boost library dir: ${Boost_LIBRARY_DIR}")
-    message("Found Boost libraries: ${Boost_LIBRARIES}")
-    include_directories(${Boost_INCLUDE_DIRS})
-endif ()
+    if (Boost_INCLUDE_DIRS)
+        if (${Boost_VERSION} VERSION_LESS 1.60.0)
+            message(FATAL_ERROR
+                    "Found Boost library is too old; required is version 1.60.0 or newer!")
+        endif ()
+        message("Found Boost include dir: ${Boost_INCLUDE_DIR}")
+        message("Found Boost library dir: ${Boost_LIBRARY_DIR}")
+        message("Found Boost libraries: ${Boost_LIBRARIES}")
+        include_directories(${Boost_INCLUDE_DIRS})
+    endif ()
+endif()
 
 #Find Python, NumPy and PyBind11
 find_package(Python COMPONENTS Interpreter Development)

guiapp/src/components/GetNearFieldSettings.vue

@@ -50,9 +50,10 @@
         <div class="col-xs-grow col-sm">
           <div class="row justify-xs-center justify-sm-start items-center">
             <div class="q-gutter-md q-py-sm q-px-xs">
-              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldType.both" :label="nearFieldType.both" />
-              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldType.Ek" :label="nearFieldType.Ek" />
-              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldType.Hk" :label="nearFieldType.Hk" />
+              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldPlane.all" label="All" />
+              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldPlane.Ek" label="Ek" />
+              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldPlane.Hk" label="Hk" />
+              <q-radio v-model="crossSection" dense size='sm' :val="nearFieldPlane.EH" label="EH" />
             </div>
           </div>
         </div>
@@ -70,7 +71,7 @@ import {
 import { useStore } from 'src/store'
 import InputWithUnits from 'components/InputWithUnits.vue'
 import { flexRowTitleStyle } from 'components/config'
-import { nearFieldType } from 'src/store/simulation-setup/state';
+import { nearFieldPlane } from 'src/store/simulation-setup/state';
 
 export default defineComponent({
 
@@ -102,11 +103,12 @@ export default defineComponent({
 
     const plotSideResolution = computed({
       get: () => $store.state.simulationSetup.gui.nearFieldSetup.plotSideResolution,
-      set: val => $store.commit('simulationSetup/setNearFieldPlotSideResolution', val)
+      // TODO: make InputWithUnits to handle integer input, so no need to use floor() in the next line.
+      set: val => $store.commit('simulationSetup/setNearFieldPlotSideResolution', Math.floor(val))
     })
 
     return { crossSection, isShowingHelpForInputWithUnits, flexRowTitleStyle,
-    relativePlotSize, maxComputeTime, plotSideResolution, nearFieldType}
+    relativePlotSize, maxComputeTime, plotSideResolution, nearFieldPlane}
   },
 })
 </script>

guiapp/src/components/RunSimulationNearField.vue

@@ -0,0 +1,198 @@
+<template>
+  <div class="row items-baseline">
+    <div class="col-xs-12 col-sm-auto text-weight-bold text-center q-px-md q-py-sm">
+      <q-tooltip
+          v-if=" $store.state.guiRuntime.safeFromWL > $store.state.simulationSetup.gui.nearFieldSetup.atWL ||
+                 $store.state.guiRuntime.safeToWL < $store.state.simulationSetup.gui.nearFieldSetup.atWL "
+          anchor="top middle" self="center middle"
+          class="bg-amber-4 text-black shadow-4">
+        Will use materials<br> spectrum range.
+      </q-tooltip>
+        <q-btn :loading="isRunning"
+               :disable="isRunning||!isNmieLoaded"
+               color="primary"
+               no-caps
+               :label="isNmieLoaded ? 'Run simulation' : 'Loading...'"
+               @click="runNearFieldSimulation">
+          <template #loading>
+            <q-spinner-gears />
+          </template>
+        </q-btn>
+    </div>
+    <div class="col-xs-grow col-sm q-px-xs">
+      <div class="row justify-xs-center justify-sm-start items-baseline">
+
+        <div class="col-auto">
+          <q-btn
+              color="primary"
+              no-caps
+              @click="saveSpectrumSimulation"
+          >Save</q-btn>
+        </div>
+
+      </div>
+    </div>
+  </div>
+</template>
+
+<script lang="ts">
+import {
+  defineComponent,
+  computed,
+  // watch,
+    onActivated,
+    nextTick
+} from 'vue'
+import { useStore } from 'src/store'
+import { getModeName, range, rangeInt} from 'components/utils'
+import { cloneDeep } from 'lodash'
+import { saveAs } from 'file-saver'
+
+export default defineComponent({
+  name: 'RunSimulationNearField',
+
+  setup() {
+    const $store = useStore()
+
+    const isRunning = computed({
+      get: ()=> $store.state.simulationSetup.nmies.nearField.isNmieRunning,
+      set: val => {
+        val ? $store.commit('simulationSetup/markNmieNearFieldAsStarted') : $store.commit('simulationSetup/markNmieNearFieldAsFinished')
+      }
+    })
+
+    const isNmieLoaded = computed(()=>{ return $store.state.simulationSetup.nmies.nearField.instance })
+
+    const atWL = computed(
+        () => $store.state.simulationSetup.current.nearFieldSetup.atWL)
+    //-------------------------------------------------------------------------//
+    //---------------------  Main  --------------------------------------------//
+    //-------------------------------------------------------------------------//
+    function runNearFieldSimulation() {
+      if (isRunning.value) {
+        console.log('Some Nmie is already running!')
+        return
+      }
+      isRunning.value = true
+      void nextTick(()=> {
+        $store.commit('simulationSetup/copySetupFromGuiToCurrent')
+
+        const host = $store.state.simulationSetup.current.hostIndex
+
+        try {
+          if (!$store.state.simulationSetup.nmies.nearField.instance) throw 'ERROR! Scattnlay module was not loaded'
+          const nmie = $store.state.simulationSetup.nmies.nearField.instance
+          const layers = cloneDeep($store.state.simulationSetup.current.layers)
+          const nmieStartedTime = performance.now()
+
+          nmie.SetWavelength(atWL.value)
+          nmie.ClearTarget()
+          for (const layer of layers) {
+            if (layer.material.nSpline) layer.n = layer.material.nSpline.at(atWL.value)
+            if (layer.material.kSpline) layer.k = layer.material.kSpline.at(atWL.value)
+            nmie.AddTargetLayerReIm(layer.layerWidth * host, layer.n / host, layer.k / host)
+          }
+
+          nmie.SetModeNmaxAndType(-1, -1)
+
+          nmie.RunFieldCalculationCartesian(
+              $store.state.simulationSetup.current.nearFieldSetup.plotSideResolution,
+              $store.state.simulationSetup.current.nearFieldSetup.relativePlotSize,
+              $store.state.simulationSetup.current.nearFieldSetup.crossSection,
+              0, 0, 0, 0
+          )
+
+
+          const nmieTotalRunTime = (performance.now()-nmieStartedTime)/1000
+          // console.log('Total simulation time:', nmieTotalRunTime, 's')
+          $store.commit('simulationSetup/setNmieNearFieldTotalRunTime', nmieTotalRunTime)
+
+        } catch (e) {
+          console.log('Some error:', e)
+        }
+        isRunning.value = false
+      })
+    }
+
+    onActivated(()=>{
+      if (isNmieLoaded.value) runNearFieldSimulation()
+    })
+
+
+
+    return { isRunning, isNmieLoaded,
+      runNearFieldSimulation,
+      saveSpectrumSimulation(){
+        const fileHeader = '# # You can open and plot this file using Python\n' +
+            '# # (without manually removing this header, it will be skipped), see example below.\n' +
+            '# import numpy as np\n' +
+            '# from matplotlib import pyplot as plt\n' +
+            '# data = np.genfromtxt(\'scattnlay-spectra.txt\', skip_header=21, names=True, delimiter=\', \')\n' +
+            '# x = data[data.dtype.names[0]] # x-axis has units\n' +
+            '# # Possible labels for spectrum data: Qsca, Qabs, Qext,\n' +
+            '# # Qsca_E_dipole, etc. (see last comment before spectra data)\n' +
+            '# a = data[\'Qsca\']\n' +
+            '# b = data[\'Qsca_E_dipole\']\n' +
+            '# c = data[\'Qsca_H_dipole\']\n' +
+            '# \n' +
+            '# plt.figure()\n' +
+            '# plt.plot(x, a, label=\'Qsca\')\n' +
+            '# plt.plot(x, b, label=\'Qsca E dipole\')\n' +
+            '# plt.plot(x, c, label=\'Qsca H dipole\')\n' +
+            '# plt.legend()\n' +
+            '# plt.xlabel(data.dtype.names[0].replace(\'_\', \', \'))\n' +
+            '# plt.ylabel(\'Normalized cross-sections\')\n' +
+            '# plt.show()\n\n'
+        let xTitle = 'x'
+        if ( $store.state.plotRuntime.spectrumPlots.layout.xaxis ) {
+          xTitle = String($store.state.plotRuntime.spectrumPlots.layout.xaxis.title)
+        }
+
+        let columnNames = '# ' + xTitle + ', Qsca, Qabs, Qext, '
+        const mode_n = rangeInt($store.state.simulationSetup.current.numberOfModesToPlot, 1);
+        const mode_types = range(0, 1);
+        for (const n of mode_n) {
+          for (const mode_type of mode_types) {
+            const modeTypeName = mode_type == 0 ? 'E' : 'H'
+            columnNames += 'Qsca_' + modeTypeName + '_' +getModeName(n)+', '
+            columnNames += 'Qabs_' + modeTypeName + '_' +getModeName(n)+', '
+            columnNames += 'Qext_' + modeTypeName + '_' +getModeName(n)+', '
+          }
+        }
+        columnNames = columnNames.slice(0, -2)
+        columnNames += '\n'
+        let body = ''
+        const WLs = $store.state.plotRuntime.WLsInUnits
+        const Qsca = $store.state.plotRuntime.Qsca
+        const Qabs = $store.state.plotRuntime.Qabs
+        const Qext = $store.state.plotRuntime.Qext
+        const Qsca_n = $store.state.plotRuntime.Qsca_n
+        const Qabs_n = $store.state.plotRuntime.Qabs_n
+        const Qext_n = $store.state.plotRuntime.Qext_n
+        for (let i = 0; i < WLs.length; ++i) {
+          let row = WLs[i].toString() + ', '
+              + Qsca[i].toString() + ', '
+              + Qabs[i].toString() + ', '
+              + Qext[i].toString() + ', '
+          for (const n of mode_n) {
+            for (const mode_type of mode_types) {
+              row += Qsca_n[mode_type][n - 1][i].toString() + ', '
+              row += Qabs_n[mode_type][n - 1][i].toString() + ', '
+              row += Qext_n[mode_type][n - 1][i].toString() + ', '
+            }
+          }
+          row = row.slice(0, -2)
+          row += '\n'
+          body += row
+        }
+
+        const scattnlaySpectra = new Blob([fileHeader+columnNames+body],
+            {type: 'text/plain;charset=utf-8',
+              endings: 'native'}  //TODO test if newline is correctly written in Windows, MacOS
+        )
+        saveAs(scattnlaySpectra, 'scattnlay-spectra.txt');
+      }
+    }
+  },
+})
+</script>

guiapp/src/nmiejs.d.ts

@@ -21,6 +21,10 @@ export class nmie_class {
                              from_Theta: number,   to_Theta: number,
                              from_Phi: number,   to_Phi: number,
                              isIgnoreAvailableNmax: number): void;
+    RunFieldCalculationCartesian( side_points: number, relative_side_length: number, plane_selected: number,
+                                  at_x: number, at_y: number, at_z: number,
+                                  isIgnoreAvailableNmax: number): void;
+
     GetFieldEabs(): number[];
     GetQsca(): number;
     GetQabs(): number;

guiapp/src/pages/Info.vue

@@ -39,6 +39,7 @@
           <li><a href="https://de.ifmo.ru/miecalculator/">Bulk particles<q-icon name="launch" style="padding-left: 0.1rem; padding-bottom: 1rem" /></a>   by Ivan Toftul.</li>
           <li><a href="https://saviot.cnrs.fr/mie/index.en.html">Bulk<q-icon name="launch" style="padding-left: 0.1rem; padding-bottom: 1rem" /></a> and <a href="https://saviot.cnrs.fr/miecoat/index.en.html">core-shell<q-icon name="launch" style="padding-left: 0.1rem; padding-bottom: 1rem" /></a> calculators by Lucien Saviot.</li>
           <li><a href="https://omlc.org/calc/mie_calc.html">Angle distribution<q-icon name="launch" style="padding-left: 0.1rem; padding-bottom: 1rem" /></a>  by Scott Prahl.</li>
+          <li><a href="http://www.philiplaven.com/mieplot.htm">Mie theory and Debye series<q-icon name="launch" style="padding-left: 0.1rem; padding-bottom: 1rem" /></a> mostly for atmospheric optical effects, by Philip Laven (Windows OS only)</li>
         </ul>
       </q-card-section>
     </q-card>

guiapp/src/pages/Near-field.vue

@@ -6,9 +6,9 @@
     <div class="q-ma-xs"/>
     <GetNearFieldSettings/>
     <div class="q-ma-xs"/>
-<!--    <RunSimulationNearField/>-->
+    <RunSimulationNearField/>
     <div class="col-auto">
-      Input result: {{$store.state.simulationSetup.gui.nearFieldSetup.crossSection}}
+      Input result: {{$store.state.simulationSetup.gui.nearFieldSetup.plotSideResolution}}
     </div>
 
     <div class="q-ma-lg"/>
@@ -26,7 +26,7 @@ import {
 import GetWlFromPlot from 'components/GetWlFromPlot.vue'
 import GetNearFieldSettings from 'components/GetNearFieldSettings.vue'
 import RunSimulationSpectrum from 'components/RunSimulationSpectrum.vue'
-// import RunSimulationNearField from 'components/RunSimulationNearField.vue'
+import RunSimulationNearField from 'components/RunSimulationNearField.vue'
 
 // import { useStore } from 'src/store'
 
@@ -34,7 +34,7 @@ import RunSimulationSpectrum from 'components/RunSimulationSpectrum.vue'
 export default defineComponent({
   name: 'NearField',
   components: {
-    // RunSimulationNearField,
+    RunSimulationNearField,
     GetWlFromPlot, RunSimulationSpectrum,
     GetNearFieldSettings,
   },

guiapp/src/store/gui-runtime/actions.ts

@@ -51,7 +51,7 @@ const actions: ActionTree<guiRuntimeStateInterface, StateInterface> = {
     let ys1:number[] = data_columns[1]
     let ys2:number[] = data_columns[1].map(()=>0)
     if (data_columns[2]) ys2 = data_columns[2]
-    const maxVal = 350
+    const maxVal = 350 // TODO move it to config.ts
     if (xs.length > maxVal) {
       const delta = Math.floor(xs.length / maxVal);
       let tmp_xs:number[] = []

guiapp/src/store/simulation-setup/mutations.ts

@@ -1,5 +1,5 @@
 import { MutationTree } from 'vuex';
-import { simulationSetupStateInterface as sssi, simulationSetup, layer, nearFieldType } from './state';
+import { simulationSetupStateInterface as sssi, simulationSetup, layer, nearFieldPlane } from './state';
 import { cloneDeep } from 'lodash'
 import { markRaw} from 'vue'
 
@@ -23,6 +23,26 @@ const mutation: MutationTree<sssi> = {
     state.nmies.spectrum.nmieTotalRunTime = val
   },
 
+  markNmieNearFieldAsStarted (state: sssi) {
+    state.nmies.nearField.isNmieRunning = true
+  },
+  markNmieNearFieldAsFinished(state: sssi) {
+    state.nmies.nearField.isNmieRunning = false
+  },
+  setNmieNearFieldTotalRunTime(state: sssi, val:number) {
+    state.nmies.nearField.nmieTotalRunTime = val
+  },
+
+  markNmieFarFieldAsStarted (state: sssi) {
+    state.nmies.farField.isNmieRunning = true
+  },
+  markNmieFarFieldAsFinished(state: sssi) {
+    state.nmies.farField.isNmieRunning = false
+  },
+  setNmieFarFieldTotalRunTime(state: sssi, val:number) {
+    state.nmies.farField.nmieTotalRunTime = val
+  },
+
   setCurrentState (state: sssi,
                    newVal: simulationSetup) {
     state.current = cloneDeep(newVal)
@@ -48,11 +68,11 @@ const mutation: MutationTree<sssi> = {
   setPlotLabel (state: sssi, val: string) {state.gui.plotLabel = val},
   setNumberOfModesToPlot  (state: sssi, val: number) {state.gui.numberOfModesToPlot  = val},
 
-  setNearFieldWL                 (state: sssi, val: number)        {state.gui.nearFieldSetup.atWL               = val},
-  setNearFieldRelativePlotSize   (state: sssi, val: number)        {state.gui.nearFieldSetup.relativePlotSize   = val},
-  setNearFieldPlotSideResolution (state: sssi, val: number)        {state.gui.nearFieldSetup.plotSideResolution = val},
-  setNearFieldCrossSection       (state: sssi, val: nearFieldType) {state.gui.nearFieldSetup.crossSection       = val},
-  setNearFieldMaxComputeTime     (state: sssi, val: number)        {state.gui.nearFieldSetup.maxComputeTime     = val},
+  setNearFieldWL                 (state: sssi, val: number)         {state.gui.nearFieldSetup.atWL               = val},
+  setNearFieldRelativePlotSize   (state: sssi, val: number)         {state.gui.nearFieldSetup.relativePlotSize   = val},
+  setNearFieldPlotSideResolution (state: sssi, val: number)         {state.gui.nearFieldSetup.plotSideResolution = val},
+  setNearFieldCrossSection       (state: sssi, val: nearFieldPlane) {state.gui.nearFieldSetup.crossSection       = val},
+  setNearFieldMaxComputeTime     (state: sssi, val: number)         {state.gui.nearFieldSetup.maxComputeTime     = val},
 
   setFarFieldWL    (state: sssi, val: number) { state.gui.farFieldWL  = val},
 

guiapp/src/store/simulation-setup/state.ts

@@ -19,17 +19,18 @@ export interface layer {
   k: number
 }
 
-export enum nearFieldType {
-  Ek = 'Ek',
-  Hk = 'Hk',
-  both = 'both'
+export enum nearFieldPlane {
+  all = -1,
+  Ek = 0,
+  Hk,
+  EH = 2
 }
 
 export interface nearFieldSetup {
   atWL:number
   relativePlotSize: number
   plotSideResolution: number
-  crossSection: nearFieldType
+  crossSection: nearFieldPlane
   maxComputeTime: number //in seconds
 }
 
@@ -74,7 +75,7 @@ function setupFactory(hostIndex = 1,
                         atWL: 619,
                         relativePlotSize: 2,
                         plotSideResolution: 64,
-                        crossSection: nearFieldType.both,
+                        crossSection: nearFieldPlane.all,
                         maxComputeTime: 5 //in seconds
                       },
 

scattnlay/main.py

@@ -31,14 +31,16 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 import numpy as np
 
-from scattnlay_mp import mie_mp as mie_mp_
-from scattnlay_dp import mie_dp
+mie_mp = None
+try:
+    from scattnlay_mp import mie_mp as mie_mp_
+    mie_mp = mie_mp_()
 
+from scattnlay_dp import mie_dp
 mie = mie_dp()
-mie_mp = mie_mp_()
 
 def scattcoeffs_(x, m, nmax=-1, pl=-1, mp=False):
-    if mp:
+    if mp and mie_mp:
         from scattnlay_mp import mie_mp as mie_
     else:
         from scattnlay_dp import mie_dp as mie_
@@ -114,7 +116,7 @@ def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
 
 
 def expancoeffs_(x, m, nmax=-1, pl=-1, mp=False):
-    if mp:
+    if mp and mie_mp:
         from scattnlay_mp import mie_mp as mie_
     else:
         from scattnlay_dp import mie_dp as mie_
@@ -199,7 +201,7 @@ def expancoeffs(x, m, nmax=-1, pl=-1, mp=False):
 
 
 def scattnlay_(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
-    if mp:
+    if mp and mie_mp:
         from scattnlay_mp import mie_mp as mie_
     else:
         from scattnlay_dp import mie_dp as mie_
@@ -284,7 +286,7 @@ def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
 
 
 def fieldnlay_(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
-    if mp:
+    if mp and mie_mp:
         from scattnlay_mp import mie_mp as mie_
     else:
         from scattnlay_dp import mie_dp as mie_

src/CMakeLists.txt

@@ -10,7 +10,11 @@ set(_scattnlay_python_sources
 
 # Define python extension
 add_library(python3-scattnlay SHARED ${_scattnlay_python_sources})
-target_link_libraries(python3-scattnlay PRIVATE Boost::headers ${Python_LIBRARIES})
+if (${ENABLE_MP} AND ${Boost_FOUND})
+    target_link_libraries(python3-scattnlay PRIVATE Boost::headers ${Python_LIBRARIES})
+else()
+    target_link_libraries(python3-scattnlay PRIVATE ${Python_LIBRARIES})
+endif()
 target_include_directories(python3-scattnlay PRIVATE ${NUMPY_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIR})
 set_target_properties(python3-scattnlay PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 
@@ -31,7 +35,9 @@ set(_scattnlay_farfield_sources
 
 # Define exe for far field calculation
 add_executable(farfield ${_scattnlay_farfield_sources})
-target_link_libraries(farfield PRIVATE Boost::headers)
+if (${ENABLE_MP} AND ${Boost_FOUND})
+    target_link_libraries(farfield PRIVATE Boost::headers)
+endif()
 set_target_properties(farfield PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 
 # Sources for near field calculation
@@ -42,11 +48,13 @@ set(_scattnlay_nearfield_sources
 
 # Define exe for near field calculation
 add_executable(nearfield ${_scattnlay_nearfield_sources})
-target_link_libraries(nearfield PRIVATE Boost::headers)
+if (${ENABLE_MP} AND ${Boost_FOUND})
+    target_link_libraries(nearfield PRIVATE Boost::headers)
+endif()
 set_target_properties(nearfield PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 
 # Rename files to match precision
-if (${ENABLE_MP})
+if (${ENABLE_MP} AND ${Boost_FOUND})
     set_property(TARGET python3-scattnlay APPEND_STRING PROPERTY OUTPUT_NAME "_mp")
     set_property(TARGET farfield APPEND_STRING PROPERTY OUTPUT_NAME "farfield-mp")
     set_property(TARGET nearfield APPEND_STRING PROPERTY OUTPUT_NAME "nearfield-mp")

src/nmie-applied-impl.hpp

@@ -278,114 +278,6 @@ namespace nmie {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  //**********************************************************************************//
-  // Function CONFRA ported from MIEV0.f (Wiscombe,1979)
-  // Ref. to NCAR Technical Notes, Wiscombe, 1979
-  /*
-c         Compute Bessel function ratio A-sub-N from its
-c         continued fraction using Lentz method
-
-c         ZINV = Reciprocal of argument of A
-
-
-c    I N T E R N A L    V A R I A B L E S
-c    ------------------------------------
-
-c    CAK      Term in continued fraction expansion of A (Eq. R25)
-c     a_k
-
-c    CAPT     Factor used in Lentz iteration for A (Eq. R27)
-c     T_k
-
-c    CNUMER   Numerator   in capT  (Eq. R28A)
-c     N_k
-c    CDENOM   Denominator in capT  (Eq. R28B)
-c     D_k
-
-c    CDTD     Product of two successive denominators of capT factors
-c                 (Eq. R34C)
-c     xi_1
-
-c    CNTN     Product of two successive numerators of capT factors
-c                 (Eq. R34B)
-c     xi_2
-
-c    EPS1     Ill-conditioning criterion
-c    EPS2     Convergence criterion
-
-c    KK       Subscript k of cAk  (Eq. R25B)
-c     k
-
-c    KOUNT    Iteration counter (used to prevent infinite looping)
-
-c    MAXIT    Max. allowed no. of iterations
-
-c    MM + 1  and - 1, alternately
-*/
-  template <typename FloatType>
-  std::complex<FloatType> MultiLayerMieApplied<FloatType>::calcD1confra(const int N, const std::complex<FloatType> z) {
-  // NTMR -> nmax_ - 1  \\TODO nmax_ ?
-    //int N = nmax_ - 1;
-    int KK, KOUNT, MAXIT = 10000, MM;
-    //    FloatType EPS1=1.0e-2;
-    FloatType EPS2=1.0e-8;
-    std::complex<FloatType> CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER;
-    std::complex<FloatType> one = std::complex<FloatType>(1.0,0.0);
-    std::complex<FloatType> ZINV = one/z;
-// c                                 ** Eq. R25a
-    std::complex<FloatType> CONFRA = static_cast<std::complex<FloatType> >(N + 1)*ZINV;   //debug ZINV
-    MM = - 1;
-    KK = 2*N +3; //debug 3
-// c                                 ** Eq. R25b, k=2
-    CAK    = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; //debug -3 ZINV
-    CDENOM = CAK;
-    CNUMER = CDENOM + one/CONFRA; //-3zinv+z
-    KOUNT  = 1;
-    //10 CONTINUE
-    do {      ++KOUNT;
-      if (KOUNT > MAXIT) {
-        printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
-        throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
-      }
-      MM *= - 1;      KK += 2;  //debug  mm=1 kk=5
-      CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; //    ** Eq. R25b //debug 5zinv
-     //  //c ** Eq. R32    Ill-conditioned case -- stride two terms instead of one
-     //  if (std::abs(CNUMER/CAK) >= EPS1 ||  std::abs(CDENOM/CAK) >= EPS1) {
-     //         //c                       ** Eq. R34
-     //         CNTN   = CAK*CNUMER + 1.0;
-     //         CDTD   = CAK*CDENOM + 1.0;
-     //         CONFRA = (CNTN/CDTD)*CONFRA; // ** Eq. R33
-     //         MM  *= - 1;        KK  += 2;
-     //         CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; // ** Eq. R25b
-     //         //c                        ** Eq. R35
-     //         CNUMER = CAK + CNUMER/CNTN;
-     //         CDENOM = CAK + CDENOM/CDTD;
-     //         ++KOUNT;
-     //         //GO TO  10
-     //         continue;
-     // } else { //c                           *** Well-conditioned case
-      {
-        CAPT   = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
-        // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
-       CONFRA = CAPT*CONFRA; // ** Eq. R26
-       //if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
-       //c                                  ** Check for convergence; Eq. R31
-       if (std::abs(CAPT.real() - 1.0) >= EPS2 ||  std::abs(CAPT.imag()) >= EPS2) {
-//c                                        ** Eq. R30
-         CNUMER = CAK + one/CNUMER;
-         CDENOM = CAK + one/CDENOM;
-         continue;
-         //GO TO  10
-       }  // end of if < eps2
-      }
-      break;
-    } while(1);
-    //if (N == 0)  printf(" return confra for z=(%g,%g)\n", ZINV.real(), ZINV.imag());
-    return CONFRA;
-  }
-  // ********************************************************************** //
-  // ********************************************************************** //
-  // ********************************************************************** //
   template <typename FloatType>
   void MultiLayerMieApplied<FloatType>::ConvertToSP() {
     this->MarkUncalculated();
@@ -418,6 +310,20 @@ c    MM + 1  and - 1, alternately
                                                                isIgnoreAvailableNmax == 0 ? false : true);
   }
 
+template <typename FloatType>
+void MultiLayerMieApplied<FloatType>::RunFieldCalculationCartesian(const int side_points,
+                                                            const double relative_side_length,
+                                                            const int plane_selected,
+                                                            const double at_x, const double at_y,
+                                                            const double at_z,
+                                                            const int isIgnoreAvailableNmax) {
+//  std::cout<<'test'<<std::endl;
+  this->MultiLayerMie<FloatType>::RunFieldCalculationCartesian( side_points, relative_side_length, plane_selected,
+                                                                at_x, at_y, at_z,
+                                                                isIgnoreAvailableNmax == 0 ? false : true);
+
+}
+
 //from https://toughengineer.github.io/demo/dsp/fft-perf/
 template <typename FloatType=double>
 emscripten::val toJSFloat64Array(const std::vector<double> &v) {

src/nmie-applied.hpp

@@ -64,6 +64,12 @@ namespace nmie {
                                   const double from_Theta, const double to_Theta,
                                   const double from_Phi, const double to_Phi,
                                   const int isIgnoreAvailableNmax);
+    void RunFieldCalculationCartesian(const int side_points,
+                                      const double relative_side_length,
+                                      const int plane_selected,
+                                      const double at_x, const double at_y,
+                                      const double at_z,
+                                      const int isIgnoreAvailableNmax);
     emscripten::val GetFieldEabs();
 
     void GetFailed();
@@ -175,7 +181,6 @@ namespace nmie {
 	            std::vector<std::complex<FloatType> >& h1np);
     void sphericalBessel(std::complex<FloatType> z, std::vector<std::complex<FloatType> >& bj,
 			             std::vector<std::complex<FloatType> >& by, std::vector<std::complex<FloatType> >& bd);
-    std::complex<FloatType> calcD1confra(int N, const std::complex<FloatType> z);
 
     FloatType wavelength_ = 1.0;
     FloatType total_radius_ = 0.0;

src/nmie-js-wrapper.cc

@@ -40,20 +40,21 @@ using namespace emscripten;
 nmie::MultiLayerMieApplied<double> ml_mie;
 
 EMSCRIPTEN_BINDINGS (c) {
-        class_<nmie::MultiLayerMieApplied<double>>("nmie")
-                .constructor<>()
-                .function("SetWavelength", &nmie::MultiLayerMieApplied<double>::SetWavelength)
-                .function("AddTargetLayerReIm",&nmie::MultiLayerMieApplied<double>::AddTargetLayerReIm)
-                .function("SetModeNmaxAndType",&nmie::MultiLayerMieApplied<double>::SetModeNmaxAndType)
-                .function("ClearTarget",&nmie::MultiLayerMieApplied<double>::ClearTarget)
-                .function("RunMieCalculation",&nmie::MultiLayerMieApplied<double>::RunMieCalculation)
-                .function("RunFieldCalculationPolar",&nmie::MultiLayerMieApplied<double>::RunFieldCalculationPolar)
-                .function("GetFieldEabs",&nmie::MultiLayerMieApplied<double>::GetFieldEabs)
-                .function("GetQsca",&nmie::MultiLayerMieApplied<double>::GetQsca)
-                .function("GetQext",&nmie::MultiLayerMieApplied<double>::GetQext)
-                .function("GetQabs",&nmie::MultiLayerMieApplied<double>::GetQabs)
-//                .function("bf",&nmie::MultiLayerMieApplied<double>::bf)
-        ;
+    class_<nmie::MultiLayerMieApplied<double>>("nmie")
+        .constructor<>()
+        .function("SetWavelength", &nmie::MultiLayerMieApplied<double>::SetWavelength)
+        .function("AddTargetLayerReIm",&nmie::MultiLayerMieApplied<double>::AddTargetLayerReIm)
+        .function("SetModeNmaxAndType",&nmie::MultiLayerMieApplied<double>::SetModeNmaxAndType)
+        .function("ClearTarget",&nmie::MultiLayerMieApplied<double>::ClearTarget)
+        .function("RunMieCalculation",&nmie::MultiLayerMieApplied<double>::RunMieCalculation)
+        .function("RunFieldCalculationPolar",&nmie::MultiLayerMieApplied<double>::RunFieldCalculationPolar)
+        .function("RunFieldCalculationCartesian",&nmie::MultiLayerMieApplied<double>::RunFieldCalculationCartesian)
+        .function("GetFieldEabs",&nmie::MultiLayerMieApplied<double>::GetFieldEabs)
+        .function("GetQsca",&nmie::MultiLayerMieApplied<double>::GetQsca)
+        .function("GetQext",&nmie::MultiLayerMieApplied<double>::GetQext)
+        .function("GetQabs",&nmie::MultiLayerMieApplied<double>::GetQabs)
+//              .function("bf",&nmie::MultiLayerMieApplied<double>::bf)
+    ;
 }
 
 //namespace nmie {

src/nmie-nearfield.hpp

@@ -454,6 +454,7 @@ namespace nmie {
     convertFieldsFromSphericalToCartesian();
   }  //  end of MultiLayerMie::RunFieldCalculation()
 
+// TODO do we really need this eval_delta()?
 template <typename FloatType>
 double eval_delta(const unsigned int steps, const double from_value, const double to_value) {
   auto delta = std::abs(from_value - to_value);
@@ -630,7 +631,39 @@ bool MultiLayerMie<FloatType>::GetFieldConvergence () {
   return convergence;
 }
 
+template <typename FloatType>
+void MultiLayerMie<FloatType>::RunFieldCalculationCartesian(const int side_points,
+                                                            const double relative_side_length,
+                                                            const int plane_selected,
+                                                            const double at_x, const double at_y,
+                                                            const double at_z,
+                                                            const bool isMarkUnconverged,
+                                                            const int nmax_in) {
+  std::cout << "params:"<<side_points<<' '<<relative_side_length<<' '<<plane_selected<<' '<<
+  at_x<<' '<<at_y<<' '<<at_z<<' '<<isMarkUnconverged<<' '<<nmax_in<<std::endl;
+//  const int planes_number = plane_selected == -1 ? 3 : 1;
+//  const int total_size = side_points*side_points*planes_number;
+//  std::vector<double> Xp(total_size), Yp(total_size), Zp(total_size);
+//  if (size_param_.size()<1) throw "Expect size_param_ to have at least one element before running a simulation";
+//  const double max_coord_value = size_param_.back() * relative_side_length;
+//  double xi = -max_coord_value, yi = -max_coord_value, zi = -max_coord_value;
+//  const double dx = max_coord_value*2/side_points;
+//  const double dy = dx, dz = dx;
+//  for (int plane = 0; plane < planes_number; ++plane) {
+//    int nx = side_points, ny = side_points, nz = side_points;
+//
+//    for (int i = 0; i < nx; i++) {
+//      for (int j = 0; j < ny; j++) {
+//        for (int k = 0; k < nz; k++) {
+//          Xp[i * ny * nz + j * nz + k] = xi + (double) i * dx;
+//          Yp[i * ny * nz + j * nz + k] = yi + (double) j * dy;
+//          Zp[i * ny * nz + j * nz + k] = zi + (double) k * dz;
+//        }
+//      }
+//    }
+//  }
 
+}
 
 }  // end of namespace nmie
 #endif  // SRC_NMIE_NEARFIELD_HPP_

src/nmie.hpp

@@ -39,9 +39,9 @@
 #include <vector>
 
 #include "nmie-precision.hpp"
-#ifdef MULTI_PRECISION
-#include <boost/math/constants/constants.hpp>
-#endif
+//#ifdef MULTI_PRECISION
+//#include <boost/math/constants/constants.hpp>
+//#endif
 
 namespace nmie {
   int ScattCoeffs(const unsigned int L, const int pl,
@@ -135,6 +135,7 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
 
   // constants for per mode evaluation
   enum Modes {kAll = -1, kElectric = 0, kMagnetic = 1};
+  enum Planes {kAllPlanes = -1, kEk = 0, kHk = 1, kEH = 2};
 
   template <typename FloatType = double>
   class MultiLayerMie {
@@ -161,13 +162,20 @@ inline std::complex<T> my_exp(const std::complex<T> &x) {
     void RunMieCalculation();
     void RunFieldCalculation(bool isMarkUnconverged=true);
 
-      void RunFieldCalculationPolar(const int outer_arc_points = 1,
+    void RunFieldCalculationPolar(const int outer_arc_points = 1,
                                   const int radius_points=1,
                                   const double from_Rho=0, const double to_Rho=static_cast<double>(1.),
                                   const double from_Theta=0, const double to_Theta=static_cast<double>(3.14159265358979323),
                                   const double from_Phi=0, const double to_Phi=static_cast<double>(3.14159265358979323),
                                   const bool isMarkUnconverged = true,
                                   int nmax_in = -1);
+    void RunFieldCalculationCartesian(const int side_points = 2,
+                                      const double relative_side_length = 2,
+                                      const int plane_selected = Planes::kAllPlanes,
+                                      const double at_x = 0, const double at_y = 0,
+                                      const double at_z = 0,
+                                      const bool isMarkUnconverged = true,
+                                      const int nmax_in = -1);
 
     void calcScattCoeffs();
     void calcExpanCoeffs();