Merge branch 'master' of https://github.com/ovidiopr/scattnlay.git

README.md

@@ -1,7 +1,8 @@
 ![output example](/doc/OutputExample.png)
 **Output example:** Field distribution inside layered Si\Ag\Si sphere
 and Poynting vector distribution in Ag sphere with poweflow lines
-calculated with Scattnlay.
+calculated with Scattnlay (scripts  field-SiAgSi-flow.py and
+field-Ag-flow.py from example section as [revision](https://github.com/ovidiopr/scattnlay/commit/57c7261705a5776f78420c1f486e929517d5f584) ).
 
 How to use scattnlay
 ====================
@@ -101,7 +102,12 @@ Papers
 2. "Reduction of scattering using thin all-dielectric shells designed by stochastic optimizer"
    Konstantin Ladutenko, Ovidio Peña-Rodríguez, Irina Melchakova, Ilya
    Yagupov, and Pavel Belov  J. Appl. Phys., vol. 116, pp. 184508,
-   2014 http://dx.doi.org/10.1063/1.4900529 
+   2014 http://dx.doi.org/10.1063/1.4900529
+
+3. "Superabsorption of light by nanoparticles" Konstantin Ladutenko,
+   Pavel Belov, Ovidio Peña-Rodríguez, Ali Mirzaei, Andrey
+   E. Miroshnichenko and Ilya V. Shadrivov  Nanoscale, 2015,7,
+   18897-18901 http://dx.doi.org/10.1039/C5NR05468K
 
 Acknowledgment
 --------------

examples/Ag-int.txt

@@ -0,0 +1,966 @@
+302.402	1.46011	2.6386
+302.803	1.46212	2.58831
+303.203	1.46401	2.53818
+303.604	1.46568	2.48835
+304.004	1.467	2.43897
+304.404	1.46787	2.3902
+304.805	1.46817	2.34218
+305.205	1.46779	2.29506
+305.606	1.46662	2.249
+306.006	1.46453	2.20414
+306.406	1.46141	2.16057
+306.807	1.45696	2.11779
+307.207	1.45083	2.07502
+307.608	1.44267	2.03148
+308.008	1.43212	1.98639
+308.408	1.41881	1.93897
+308.809	1.4024	1.88843
+309.209	1.38252	1.834
+309.61	1.35882	1.77489
+310.01	1.33095	1.71031
+310.41	1.29904	1.64055
+310.811	1.26512	1.56977
+311.211	1.23162	1.50305
+311.612	1.20098	1.44543
+312.012	1.174	1.39916
+312.412	1.14862	1.36144
+312.813	1.12246	1.32896
+313.213	1.09313	1.2984
+313.614	1.05826	1.26645
+314.014	1.0155	1.22981
+314.414	0.964964	1.18761
+314.815	0.911161	1.14324
+315.215	0.859034	1.1005
+315.616	0.813496	1.0632
+316.016	0.776525	1.0332
+316.416	0.744974	1.00901
+316.817	0.715181	0.98879
+317.217	0.683483	0.970693
+317.618	0.646216	0.952884
+318.018	0.59976	0.933539
+318.418	0.543756	0.912113
+318.819	0.483461	0.89026
+319.219	0.424687	0.869849
+319.62	0.373206	0.852736
+320.02	0.331754	0.839738
+320.42	0.297907	0.829895
+320.821	0.268739	0.822078
+321.221	0.241324	0.815156
+321.622	0.212735	0.807999
+322.022	0.180046	0.799477
+322.422	0.141494	0.788869
+322.823	0.0992597	0.776848
+323.223	0.0563613	0.764383
+323.624	0.0158164	0.752439
+324.024	-0.0197298	0.741868
+324.424	-0.0503055	0.73268
+324.825	-0.0771062	0.724522
+325.225	-0.101332	0.717034
+325.626	-0.124185	0.709861
+326.026	-0.146863	0.702646
+326.426	-0.170562	0.695033
+326.827	-0.196039	0.686842
+327.227	-0.223341	0.678171
+327.628	-0.25245	0.669143
+328.028	-0.283348	0.659884
+328.428	-0.316018	0.650517
+328.829	-0.350442	0.641167
+329.229	-0.386603	0.631959
+329.63	-0.424481	0.623015
+330.03	-0.46406	0.614461
+330.43	-0.505323	0.60642
+330.831	-0.548204	0.599019
+331.231	-0.591893	0.592402
+331.632	-0.634963	0.586733
+332.032	-0.675971	0.582173
+332.432	-0.713473	0.578885
+332.833	-0.746512	0.576943
+333.233	-0.775661	0.57614
+333.634	-0.801793	0.576217
+334.034	-0.825781	0.576912
+334.434	-0.848497	0.577964
+334.835	-0.870814	0.579113
+335.235	-0.893598	0.5801
+335.636	-0.917345	0.580781
+336.036	-0.941978	0.581187
+336.436	-0.967374	0.581366
+336.837	-0.993411	0.581363
+337.237	-1.01997	0.581226
+337.638	-1.04691	0.581
+338.038	-1.07413	0.580733
+338.438	-1.1015	0.580471
+338.839	-1.12889	0.58026
+339.239	-1.15618	0.580148
+339.64	-1.18325	0.58018
+340.04	-1.20998	0.580396
+340.44	-1.23637	0.580794
+340.841	-1.26241	0.581353
+341.241	-1.28812	0.582053
+341.642	-1.31353	0.582874
+342.042	-1.33862	0.583796
+342.442	-1.36343	0.584799
+342.843	-1.38797	0.585862
+343.243	-1.41224	0.586966
+343.644	-1.43626	0.588091
+344.044	-1.46004	0.589216
+344.444	-1.4836	0.590321
+344.845	-1.50695	0.591391
+345.245	-1.53011	0.592422
+345.646	-1.55308	0.593411
+346.046	-1.57589	0.594359
+346.446	-1.59855	0.595263
+346.847	-1.62107	0.596121
+347.247	-1.64348	0.596932
+347.648	-1.66579	0.597694
+348.048	-1.68801	0.598406
+348.448	-1.71016	0.599067
+348.849	-1.73226	0.599674
+349.249	-1.75431	0.600226
+349.65	-1.77634	0.600722
+350.05	-1.79836	0.60116
+350.45	-1.82039	0.601538
+350.851	-1.84245	0.601855
+351.251	-1.86454	0.602109
+351.652	-1.88668	0.602299
+352.052	-1.90889	0.602423
+352.452	-1.93119	0.60248
+352.853	-1.95359	0.602468
+353.253	-1.97611	0.602385
+353.654	-1.99876	0.60223
+354.054	-2.02155	0.602002
+354.454	-2.04451	0.601699
+354.855	-2.0676	0.601326
+355.255	-2.0908	0.600895
+355.656	-2.11404	0.600414
+356.056	-2.1373	0.599895
+356.456	-2.16052	0.599348
+356.857	-2.18366	0.598783
+357.257	-2.20668	0.598209
+357.658	-2.22953	0.597638
+358.058	-2.25218	0.59708
+358.458	-2.27456	0.596544
+358.859	-2.29666	0.596041
+359.259	-2.31841	0.595582
+359.66	-2.33977	0.595176
+360.06	-2.36071	0.594834
+360.46	-2.38117	0.594566
+360.861	-2.40111	0.594382
+361.261	-2.4205	0.594293
+361.662	-2.43928	0.594308
+362.062	-2.45741	0.594438
+362.462	-2.47484	0.594693
+362.863	-2.49155	0.595084
+363.263	-2.50747	0.59562
+363.664	-2.52256	0.596312
+364.064	-2.53679	0.59717
+364.464	-2.55011	0.598205
+364.865	-2.56247	0.599426
+365.265	-2.57388	0.600834
+365.666	-2.58441	0.602419
+366.066	-2.5941	0.604171
+366.466	-2.60304	0.606078
+366.867	-2.61128	0.60813
+367.267	-2.6189	0.610316
+367.668	-2.62594	0.612624
+368.068	-2.63249	0.615045
+368.468	-2.6386	0.617566
+368.869	-2.64435	0.620178
+369.269	-2.64978	0.622869
+369.67	-2.65498	0.625629
+370.07	-2.66	0.628446
+370.47	-2.66491	0.63131
+370.871	-2.66978	0.63421
+371.271	-2.67467	0.637135
+371.672	-2.67964	0.640075
+372.072	-2.68476	0.643017
+372.472	-2.6901	0.645952
+372.873	-2.69571	0.648868
+373.273	-2.70167	0.651755
+373.674	-2.70804	0.654602
+374.074	-2.71488	0.657398
+374.474	-2.72227	0.660132
+374.875	-2.73026	0.662793
+375.275	-2.73891	0.66537
+375.676	-2.7483	0.667852
+376.076	-2.75848	0.670231
+376.476	-2.76943	0.672504
+376.877	-2.78111	0.674672
+377.277	-2.7935	0.676735
+377.678	-2.80655	0.678693
+378.078	-2.82023	0.680546
+378.478	-2.83451	0.682296
+378.879	-2.84936	0.683943
+379.279	-2.86473	0.685486
+379.68	-2.88059	0.686926
+380.08	-2.89691	0.688265
+380.48	-2.91365	0.689501
+380.881	-2.93079	0.690635
+381.281	-2.94827	0.691669
+381.682	-2.96608	0.692601
+382.082	-2.98417	0.693433
+382.482	-3.00251	0.694165
+382.883	-3.02106	0.694797
+383.283	-3.03979	0.695329
+383.684	-3.05867	0.695763
+384.084	-3.07766	0.696098
+384.484	-3.09673	0.696335
+384.885	-3.11583	0.696474
+385.285	-3.13495	0.696515
+385.686	-3.15403	0.696459
+386.086	-3.17305	0.696307
+386.486	-3.19197	0.696058
+386.887	-3.21076	0.695713
+387.287	-3.22939	0.695273
+387.688	-3.24781	0.694737
+388.088	-3.26602	0.694111
+388.488	-3.28402	0.693402
+388.889	-3.30184	0.692618
+389.289	-3.31948	0.691769
+389.69	-3.33696	0.690863
+390.09	-3.35429	0.689908
+390.49	-3.37149	0.688912
+390.891	-3.38856	0.687885
+391.291	-3.40552	0.686833
+391.692	-3.42238	0.685767
+392.092	-3.43916	0.684695
+392.492	-3.45587	0.683624
+392.893	-3.47252	0.682563
+393.293	-3.48912	0.681521
+393.694	-3.5057	0.680506
+394.094	-3.52225	0.679527
+394.494	-3.53879	0.678592
+394.895	-3.55535	0.677709
+395.295	-3.57192	0.676888
+395.696	-3.58852	0.676136
+396.096	-3.60518	0.675461
+396.496	-3.62189	0.674873
+396.897	-3.63867	0.67438
+397.297	-3.65554	0.67399
+397.698	-3.6725	0.673711
+398.098	-3.68958	0.673553
+398.498	-3.70678	0.673523
+398.899	-3.72411	0.67363
+399.299	-3.7416	0.673883
+399.7	-3.75925	0.674289
+400.1	-3.77707	0.674858
+400.501	-3.79508	0.675592
+400.901	-3.81327	0.676482
+401.301	-3.83162	0.67752
+401.702	-3.85013	0.678696
+402.102	-3.8688	0.68
+402.503	-3.88761	0.681423
+402.903	-3.90656	0.682955
+403.303	-3.92565	0.684587
+403.704	-3.94485	0.68631
+404.104	-3.96417	0.688114
+404.505	-3.9836	0.689989
+404.905	-4.00313	0.691926
+405.305	-4.02275	0.693915
+405.706	-4.04245	0.695948
+406.106	-4.06224	0.698014
+406.507	-4.08209	0.700104
+406.907	-4.102	0.702209
+407.307	-4.12197	0.704319
+407.708	-4.14199	0.706425
+408.108	-4.16204	0.708516
+408.509	-4.18213	0.710585
+408.909	-4.20224	0.71262
+409.309	-4.22237	0.714613
+409.71	-4.2425	0.716555
+410.11	-4.26264	0.718435
+410.511	-4.28277	0.720244
+410.911	-4.30288	0.721974
+411.311	-4.32298	0.723613
+411.712	-4.34304	0.725153
+412.112	-4.36307	0.726585
+412.513	-4.38306	0.727898
+412.913	-4.40299	0.729084
+413.313	-4.42286	0.730133
+413.714	-4.44267	0.731038
+414.114	-4.4624	0.731803
+414.515	-4.48208	0.732436
+414.915	-4.50168	0.732944
+415.315	-4.52121	0.733332
+415.716	-4.54068	0.733608
+416.116	-4.56008	0.733779
+416.517	-4.57941	0.733851
+416.917	-4.59867	0.733831
+417.317	-4.61786	0.733726
+417.718	-4.63698	0.733542
+418.118	-4.65604	0.733287
+418.519	-4.67502	0.732967
+418.919	-4.69393	0.732588
+419.319	-4.71277	0.732158
+419.72	-4.73155	0.731684
+420.12	-4.75025	0.731172
+420.521	-4.76888	0.730628
+420.921	-4.78744	0.730061
+421.321	-4.80592	0.729476
+421.722	-4.82434	0.72888
+422.122	-4.84269	0.72828
+422.523	-4.86096	0.727682
+422.923	-4.87916	0.727094
+423.323	-4.89729	0.726523
+423.724	-4.91534	0.725974
+424.124	-4.93332	0.725455
+424.525	-4.95123	0.724973
+424.925	-4.96907	0.724534
+425.325	-4.98683	0.724145
+425.726	-5.00452	0.723813
+426.126	-5.02213	0.723545
+426.527	-5.03967	0.723347
+426.927	-5.05713	0.723226
+427.327	-5.07452	0.723188
+427.728	-5.09184	0.723242
+428.128	-5.10908	0.723386
+428.529	-5.12626	0.723617
+428.929	-5.14337	0.72393
+429.329	-5.16043	0.724322
+429.73	-5.17743	0.724788
+430.13	-5.19439	0.725324
+430.531	-5.2113	0.725926
+430.931	-5.22818	0.726588
+431.331	-5.24503	0.727308
+431.732	-5.26184	0.72808
+432.132	-5.27864	0.728901
+432.533	-5.29542	0.729766
+432.933	-5.31218	0.730671
+433.333	-5.32894	0.731611
+433.734	-5.34569	0.732583
+434.134	-5.36245	0.733582
+434.535	-5.37922	0.734604
+434.935	-5.396	0.735644
+435.335	-5.41279	0.736699
+435.736	-5.42961	0.737763
+436.136	-5.44645	0.738834
+436.537	-5.46333	0.739906
+436.937	-5.48025	0.740975
+437.337	-5.49721	0.742037
+437.738	-5.51421	0.743088
+438.138	-5.53127	0.744123
+438.539	-5.54838	0.745138
+438.939	-5.56556	0.746129
+439.339	-5.58281	0.747092
+439.74	-5.60012	0.748022
+440.14	-5.61752	0.748916
+440.541	-5.635	0.749768
+440.941	-5.65256	0.750574
+441.341	-5.67022	0.751331
+441.742	-5.68797	0.752034
+442.142	-5.70583	0.752679
+442.543	-5.72379	0.753261
+442.943	-5.74187	0.753776
+443.343	-5.76006	0.754223
+443.744	-5.77836	0.754602
+444.144	-5.79676	0.754917
+444.545	-5.81527	0.755168
+444.945	-5.83388	0.755358
+445.345	-5.85259	0.755488
+445.746	-5.87139	0.755561
+446.146	-5.89028	0.755577
+446.547	-5.90926	0.75554
+446.947	-5.92832	0.75545
+447.347	-5.94746	0.755311
+447.748	-5.96667	0.755123
+448.148	-5.98595	0.754888
+448.549	-6.00531	0.754609
+448.949	-6.02472	0.754287
+449.349	-6.0442	0.753923
+449.75	-6.06374	0.753521
+450.15	-6.08333	0.753082
+450.551	-6.10297	0.752607
+450.951	-6.12266	0.752099
+451.351	-6.1424	0.751559
+451.752	-6.16217	0.75099
+452.152	-6.18198	0.750392
+452.553	-6.20182	0.749769
+452.953	-6.2217	0.749121
+453.353	-6.24159	0.748451
+453.754	-6.26152	0.74776
+454.154	-6.28146	0.747051
+454.555	-6.30141	0.746326
+454.955	-6.32138	0.745585
+455.355	-6.34136	0.744831
+455.756	-6.36134	0.744067
+456.156	-6.38132	0.743293
+456.557	-6.4013	0.742511
+456.957	-6.42127	0.741724
+457.357	-6.44124	0.740933
+457.758	-6.46119	0.740141
+458.158	-6.48113	0.739348
+458.559	-6.50104	0.738558
+458.959	-6.52094	0.737771
+459.359	-6.5408	0.73699
+459.76	-6.56064	0.736216
+460.16	-6.58045	0.735449
+460.561	-6.60023	0.734691
+460.961	-6.61998	0.733941
+461.361	-6.6397	0.733201
+461.762	-6.65939	0.732471
+462.162	-6.67905	0.731751
+462.563	-6.69868	0.731043
+462.963	-6.71828	0.730346
+463.363	-6.73786	0.729662
+463.764	-6.7574	0.728991
+464.164	-6.77691	0.728334
+464.565	-6.79639	0.727691
+464.965	-6.81584	0.727062
+465.365	-6.83526	0.726449
+465.766	-6.85465	0.725853
+466.166	-6.874	0.725273
+466.567	-6.89333	0.72471
+466.967	-6.91262	0.724165
+467.367	-6.93189	0.723638
+467.768	-6.95112	0.723131
+468.168	-6.97032	0.722643
+468.569	-6.98949	0.722176
+468.969	-7.00863	0.72173
+469.369	-7.02774	0.721305
+469.77	-7.04681	0.720902
+470.17	-7.06585	0.720522
+470.571	-7.08486	0.720165
+470.971	-7.10384	0.719832
+471.371	-7.12279	0.719524
+471.772	-7.1417	0.71924
+472.172	-7.16058	0.718983
+472.573	-7.17943	0.718752
+472.973	-7.19824	0.718548
+473.373	-7.21702	0.718371
+473.774	-7.23577	0.718223
+474.174	-7.25449	0.718103
+474.575	-7.27317	0.718013
+474.975	-7.29182	0.717952
+475.375	-7.31043	0.717923
+475.776	-7.32901	0.717924
+476.176	-7.34756	0.717957
+476.577	-7.36607	0.718023
+476.977	-7.38455	0.718122
+477.377	-7.403	0.718253
+477.778	-7.42141	0.718418
+478.178	-7.4398	0.718614
+478.579	-7.45816	0.718842
+478.979	-7.47651	0.719101
+479.379	-7.49484	0.719389
+479.78	-7.51315	0.719707
+480.18	-7.53145	0.720053
+480.581	-7.54975	0.720427
+480.981	-7.56804	0.720829
+481.381	-7.58633	0.721256
+481.782	-7.60463	0.72171
+482.182	-7.62293	0.722189
+482.583	-7.64123	0.722692
+482.983	-7.65956	0.723219
+483.383	-7.67789	0.72377
+483.784	-7.69625	0.724342
+484.184	-7.71463	0.724937
+484.585	-7.73303	0.725552
+484.985	-7.75147	0.726188
+485.385	-7.76993	0.726843
+485.786	-7.78843	0.727518
+486.186	-7.80698	0.728211
+486.587	-7.82556	0.728921
+486.987	-7.84419	0.729648
+487.387	-7.86286	0.730392
+487.788	-7.88159	0.731151
+488.188	-7.90037	0.731926
+488.589	-7.91922	0.732714
+488.989	-7.93812	0.733516
+489.389	-7.95709	0.734331
+489.79	-7.97613	0.735158
+490.19	-7.99524	0.735996
+490.591	-8.01442	0.736846
+490.991	-8.03369	0.737705
+491.391	-8.05303	0.738574
+491.792	-8.07246	0.739451
+492.192	-8.09198	0.740337
+492.593	-8.11159	0.74123
+492.993	-8.13129	0.74213
+493.393	-8.1511	0.743036
+493.794	-8.171	0.743947
+494.194	-8.19101	0.744862
+494.595	-8.21113	0.745782
+494.995	-8.23135	0.746705
+495.395	-8.2517	0.747631
+495.796	-8.27216	0.748558
+496.196	-8.29274	0.749487
+496.597	-8.31344	0.750417
+496.997	-8.33425	0.751348
+497.397	-8.35518	0.75228
+497.798	-8.37621	0.753213
+498.198	-8.39734	0.754147
+498.599	-8.41857	0.755082
+498.999	-8.4399	0.756019
+499.399	-8.46132	0.756956
+499.8	-8.48282	0.757894
+500.2	-8.5044	0.758834
+500.601	-8.52607	0.759774
+501.001	-8.5478	0.760716
+501.401	-8.56961	0.761658
+501.802	-8.59148	0.762602
+502.202	-8.61341	0.763546
+502.603	-8.6354	0.764492
+503.003	-8.65744	0.765439
+503.403	-8.67953	0.766387
+503.804	-8.70166	0.767335
+504.204	-8.72384	0.768285
+504.605	-8.74605	0.769236
+505.005	-8.7683	0.770188
+505.405	-8.79057	0.771141
+505.806	-8.81286	0.772095
+506.206	-8.83518	0.77305
+506.607	-8.85751	0.774006
+507.007	-8.87986	0.774963
+507.407	-8.90221	0.775921
+507.808	-8.92456	0.77688
+508.208	-8.94692	0.777841
+508.609	-8.96927	0.778802
+509.009	-8.99161	0.779764
+509.409	-9.01393	0.780727
+509.81	-9.03624	0.781692
+510.21	-9.05853	0.782657
+510.611	-9.0808	0.783623
+511.011	-9.10303	0.784591
+511.411	-9.12523	0.785559
+511.812	-9.14739	0.786528
+512.212	-9.16951	0.787499
+512.613	-9.19159	0.78847
+513.013	-9.21361	0.789443
+513.413	-9.23558	0.790416
+513.814	-9.25749	0.791391
+514.214	-9.27934	0.792366
+514.615	-9.30112	0.793343
+515.015	-9.32283	0.79432
+515.415	-9.34447	0.795299
+515.816	-9.36603	0.796279
+516.216	-9.3875	0.797259
+516.617	-9.40889	0.798241
+517.017	-9.43018	0.799223
+517.417	-9.45139	0.800206
+517.818	-9.47251	0.801189
+518.218	-9.49355	0.802171
+518.619	-9.51451	0.803152
+519.019	-9.53539	0.804131
+519.419	-9.5562	0.805108
+519.82	-9.57694	0.806082
+520.22	-9.59761	0.807052
+520.621	-9.61821	0.808018
+521.021	-9.63874	0.80898
+521.421	-9.65922	0.809936
+521.822	-9.67964	0.810887
+522.222	-9.70001	0.811832
+522.623	-9.72032	0.81277
+523.023	-9.74058	0.8137
+523.423	-9.7608	0.814623
+523.824	-9.78098	0.815537
+524.224	-9.80111	0.816441
+524.625	-9.8212	0.817337
+525.025	-9.84126	0.818222
+525.425	-9.86129	0.819097
+525.826	-9.88128	0.81996
+526.226	-9.90125	0.820811
+526.627	-9.9212	0.82165
+527.027	-9.94112	0.822476
+527.427	-9.96103	0.823289
+527.828	-9.98092	0.824087
+528.228	-10.0008	0.824871
+528.629	-10.0207	0.82564
+529.029	-10.0405	0.826393
+529.429	-10.0604	0.827129
+529.83	-10.0802	0.827849
+530.23	-10.1001	0.828552
+530.631	-10.1199	0.829236
+531.031	-10.1398	0.829902
+531.431	-10.1597	0.830549
+531.832	-10.1796	0.831176
+532.232	-10.1995	0.831783
+532.633	-10.2194	0.83237
+533.033	-10.2393	0.832935
+533.433	-10.2593	0.833478
+533.834	-10.2793	0.833999
+534.234	-10.2993	0.834497
+534.635	-10.3193	0.834971
+535.035	-10.3394	0.835421
+535.435	-10.3595	0.835847
+535.836	-10.3797	0.836247
+536.236	-10.3999	0.836621
+536.637	-10.4202	0.83697
+537.037	-10.4405	0.837291
+537.437	-10.4608	0.837585
+537.838	-10.4812	0.837851
+538.238	-10.5017	0.838088
+538.639	-10.5222	0.838297
+539.039	-10.5428	0.838475
+539.439	-10.5634	0.838624
+539.84	-10.5841	0.838743
+540.24	-10.6049	0.838834
+540.641	-10.6257	0.838897
+541.041	-10.6466	0.838933
+541.441	-10.6676	0.838944
+541.842	-10.6886	0.83893
+542.242	-10.7097	0.838892
+542.643	-10.7308	0.838832
+543.043	-10.752	0.838749
+543.443	-10.7732	0.838646
+543.844	-10.7945	0.838523
+544.244	-10.8159	0.83838
+544.645	-10.8373	0.83822
+545.045	-10.8587	0.838043
+545.445	-10.8802	0.83785
+545.846	-10.9018	0.837641
+546.246	-10.9234	0.837418
+546.647	-10.945	0.837182
+547.047	-10.9667	0.836934
+547.447	-10.9884	0.836674
+547.848	-11.0102	0.836404
+548.248	-11.032	0.836125
+548.649	-11.0539	0.835837
+549.049	-11.0757	0.835542
+549.449	-11.0977	0.83524
+549.85	-11.1196	0.834932
+550.25	-11.1416	0.83462
+550.651	-11.1637	0.834304
+551.051	-11.1857	0.833986
+551.451	-11.2078	0.833665
+551.852	-11.23	0.833344
+552.252	-11.2521	0.833023
+552.653	-11.2743	0.832703
+553.053	-11.2965	0.832385
+553.453	-11.3188	0.83207
+553.854	-11.3411	0.83176
+554.254	-11.3634	0.831453
+554.655	-11.3857	0.831153
+555.055	-11.408	0.83086
+555.455	-11.4304	0.830574
+555.856	-11.4528	0.830298
+556.256	-11.4752	0.830031
+556.657	-11.4976	0.829774
+557.057	-11.52	0.829529
+557.457	-11.5425	0.829297
+557.858	-11.565	0.829078
+558.258	-11.5874	0.828874
+558.659	-11.6099	0.828685
+559.059	-11.6324	0.828512
+559.459	-11.6549	0.828357
+559.86	-11.6775	0.82822
+560.26	-11.7	0.828103
+560.661	-11.7225	0.828005
+561.061	-11.7451	0.827929
+561.461	-11.7676	0.827875
+561.862	-11.7902	0.827844
+562.262	-11.8127	0.827837
+562.663	-11.8353	0.827855
+563.063	-11.8579	0.827899
+563.463	-11.8804	0.82797
+563.864	-11.903	0.828068
+564.264	-11.9255	0.828195
+564.665	-11.9481	0.828349
+565.065	-11.9706	0.828531
+565.465	-11.9931	0.82874
+565.866	-12.0157	0.828977
+566.266	-12.0382	0.82924
+566.667	-12.0607	0.82953
+567.067	-12.0833	0.829846
+567.467	-12.1058	0.830188
+567.868	-12.1283	0.830556
+568.268	-12.1508	0.83095
+568.669	-12.1733	0.83137
+569.069	-12.1958	0.831814
+569.469	-12.2183	0.832284
+569.87	-12.2408	0.832778
+570.27	-12.2633	0.833297
+570.671	-12.2858	0.83384
+571.071	-12.3082	0.834407
+571.471	-12.3307	0.834998
+571.872	-12.3532	0.835613
+572.272	-12.3756	0.83625
+572.673	-12.3981	0.836911
+573.073	-12.4205	0.837595
+573.473	-12.4429	0.838302
+573.874	-12.4653	0.83903
+574.274	-12.4878	0.839781
+574.675	-12.5102	0.840554
+575.075	-12.5326	0.841349
+575.475	-12.5549	0.842165
+575.876	-12.5773	0.843002
+576.276	-12.5997	0.84386
+576.677	-12.622	0.84474
+577.077	-12.6444	0.845639
+577.477	-12.6667	0.846559
+577.878	-12.6891	0.847499
+578.278	-12.7114	0.848458
+578.679	-12.7337	0.849438
+579.079	-12.756	0.850436
+579.479	-12.7783	0.851454
+579.88	-12.8005	0.852491
+580.28	-12.8228	0.853546
+580.681	-12.8451	0.85462
+581.081	-12.8673	0.855712
+581.481	-12.8895	0.856821
+581.882	-12.9117	0.857949
+582.282	-12.9339	0.859094
+582.683	-12.9561	0.860256
+583.083	-12.9783	0.861435
+583.483	-13.0005	0.862631
+583.884	-13.0226	0.863844
+584.284	-13.0448	0.865073
+584.685	-13.0669	0.866318
+585.085	-13.089	0.867578
+585.485	-13.1111	0.868855
+585.886	-13.1332	0.870146
+586.286	-13.1552	0.871453
+586.687	-13.1773	0.872775
+587.087	-13.1993	0.874111
+587.487	-13.2213	0.875462
+587.888	-13.2433	0.876827
+588.288	-13.2653	0.878206
+588.689	-13.2873	0.879599
+589.089	-13.3092	0.881005
+589.489	-13.3312	0.882424
+589.89	-13.3531	0.883857
+590.29	-13.375	0.885302
+590.691	-13.3969	0.88676
+591.091	-13.4188	0.88823
+591.491	-13.4406	0.889712
+591.892	-13.4625	0.891207
+592.292	-13.4843	0.892713
+592.693	-13.5061	0.894231
+593.093	-13.5279	0.895761
+593.493	-13.5497	0.897302
+593.894	-13.5715	0.898854
+594.294	-13.5933	0.900418
+594.695	-13.6151	0.901993
+595.095	-13.6369	0.903578
+595.495	-13.6586	0.905174
+595.896	-13.6804	0.906781
+596.296	-13.7022	0.908398
+596.697	-13.724	0.910026
+597.097	-13.7458	0.911663
+597.497	-13.7676	0.91331
+597.898	-13.7894	0.914967
+598.298	-13.8112	0.916634
+598.699	-13.833	0.91831
+599.099	-13.8548	0.919996
+599.499	-13.8767	0.921691
+599.9	-13.8985	0.923394
+600.3	-13.9204	0.925107
+600.701	-13.9423	0.926828
+601.101	-13.9642	0.928558
+601.502	-13.9861	0.930296
+601.902	-14.0081	0.932043
+602.302	-14.03	0.933798
+602.703	-14.052	0.93556
+603.103	-14.074	0.93733
+603.504	-14.0961	0.939108
+603.904	-14.1181	0.940894
+604.304	-14.1402	0.942687
+604.705	-14.1623	0.944487
+605.105	-14.1845	0.946294
+605.506	-14.2067	0.948108
+605.906	-14.2289	0.949928
+606.306	-14.2511	0.951755
+606.707	-14.2734	0.953589
+607.107	-14.2958	0.955429
+607.508	-14.3181	0.957275
+607.908	-14.3405	0.959127
+608.308	-14.363	0.960985
+608.709	-14.3855	0.962848
+609.109	-14.408	0.964717
+609.51	-14.4306	0.966591
+609.91	-14.4532	0.96847
+610.31	-14.4759	0.970355
+610.711	-14.4987	0.972244
+611.111	-14.5214	0.974138
+611.512	-14.5443	0.976037
+611.912	-14.5672	0.97794
+612.312	-14.5901	0.979848
+612.713	-14.6131	0.981759
+613.113	-14.6362	0.983675
+613.514	-14.6593	0.985594
+613.914	-14.6825	0.987517
+614.314	-14.7057	0.989443
+614.715	-14.729	0.991373
+615.115	-14.7524	0.993306
+615.516	-14.7758	0.995242
+615.916	-14.7993	0.997181
+616.316	-14.8229	0.999123
+616.717	-14.8466	1.00107
+617.117	-14.8703	1.00301
+617.518	-14.8941	1.00496
+617.918	-14.9179	1.00691
+618.318	-14.9419	1.00887
+618.719	-14.9659	1.01082
+619.119	-14.99	1.01278
+619.52	-15.0142	1.01474
+619.92	-15.0384	1.0167
+620.32	-15.0627	1.01866
+620.721	-15.0872	1.02062
+621.121	-15.1117	1.02258
+621.522	-15.1362	1.02455
+621.922	-15.1609	1.02651
+622.322	-15.1856	1.02847
+622.723	-15.2105	1.03044
+623.123	-15.2353	1.0324
+623.524	-15.2603	1.03436
+623.924	-15.2853	1.03633
+624.324	-15.3104	1.03829
+624.725	-15.3356	1.04025
+625.125	-15.3609	1.04221
+625.526	-15.3862	1.04417
+625.926	-15.4116	1.04612
+626.326	-15.437	1.04808
+626.727	-15.4625	1.05003
+627.127	-15.4881	1.05198
+627.528	-15.5138	1.05393
+627.928	-15.5395	1.05587
+628.328	-15.5652	1.05781
+628.729	-15.5911	1.05975
+629.129	-15.617	1.06169
+629.53	-15.6429	1.06362
+629.93	-15.6689	1.06555
+630.33	-15.695	1.06747
+630.731	-15.7211	1.06939
+631.131	-15.7473	1.07131
+631.532	-15.7735	1.07322
+631.932	-15.7998	1.07512
+632.332	-15.8261	1.07702
+632.733	-15.8525	1.07892
+633.133	-15.879	1.08081
+633.534	-15.9054	1.0827
+633.934	-15.932	1.08458
+634.334	-15.9585	1.08645
+634.735	-15.9852	1.08832
+635.135	-16.0118	1.09018
+635.536	-16.0385	1.09203
+635.936	-16.0653	1.09388
+636.336	-16.0921	1.09572
+636.737	-16.1189	1.09755
+637.137	-16.1458	1.09937
+637.538	-16.1727	1.10119
+637.938	-16.1996	1.103
+638.338	-16.2266	1.1048
+638.739	-16.2536	1.1066
+639.139	-16.2807	1.10838
+639.54	-16.3078	1.11016
+639.94	-16.3349	1.11192
+640.34	-16.362	1.11368
+640.741	-16.3892	1.11543
+641.141	-16.4164	1.11717
+641.542	-16.4437	1.1189
+641.942	-16.4709	1.12062
+642.342	-16.4982	1.12233
+642.743	-16.5255	1.12403
+643.143	-16.5529	1.12571
+643.544	-16.5802	1.12739
+643.944	-16.6076	1.12906
+644.344	-16.635	1.13071
+644.745	-16.6624	1.13236
+645.145	-16.6899	1.13399
+645.546	-16.7173	1.13561
+645.946	-16.7448	1.13722
+646.346	-16.7723	1.13882
+646.747	-16.7998	1.1404
+647.147	-16.8273	1.14197
+647.548	-16.8549	1.14353
+647.948	-16.8824	1.14507
+648.348	-16.91	1.14661
+648.749	-16.9375	1.14812
+649.149	-16.9651	1.14963
+649.55	-16.9927	1.15112
+649.95	-17.0202	1.1526
+650.35	-17.0478	1.15406
+650.751	-17.0754	1.15551
+651.151	-17.103	1.15694
+651.552	-17.1306	1.15836
+651.952	-17.1582	1.15976
+652.352	-17.1858	1.16115
+652.753	-17.2134	1.16252
+653.153	-17.241	1.16388
+653.554	-17.2686	1.16522
+653.954	-17.2962	1.16655
+654.354	-17.3238	1.16786
+654.755	-17.3514	1.16915
+655.155	-17.379	1.17043
+655.556	-17.4066	1.1717
+655.956	-17.4342	1.17295
+656.356	-17.4618	1.17419
+656.757	-17.4893	1.17542
+657.157	-17.5169	1.17663
+657.558	-17.5445	1.17782
+657.958	-17.5721	1.179
+658.358	-17.5997	1.18017
+658.759	-17.6273	1.18133
+659.159	-17.6548	1.18247
+659.56	-17.6824	1.1836
+659.96	-17.71	1.18472
+660.36	-17.7376	1.18583
+660.761	-17.7651	1.18692
+661.161	-17.7927	1.188
+661.562	-17.8203	1.18907
+661.962	-17.8478	1.19013
+662.362	-17.8754	1.19117
+662.763	-17.903	1.19221
+663.163	-17.9305	1.19323
+663.564	-17.9581	1.19424
+663.964	-17.9857	1.19524
+664.364	-18.0132	1.19623
+664.765	-18.0408	1.19721
+665.165	-18.0684	1.19818
+665.566	-18.0959	1.19914
+665.966	-18.1235	1.20009
+666.366	-18.151	1.20103
+666.767	-18.1786	1.20196
+667.167	-18.2061	1.20288
+667.568	-18.2337	1.20379
+667.968	-18.2612	1.20469
+668.368	-18.2888	1.20558
+668.769	-18.3164	1.20647
+669.169	-18.3439	1.20734
+669.57	-18.3715	1.20821
+669.97	-18.399	1.20906
+670.37	-18.4266	1.20991
+670.771	-18.4541	1.21075
+671.171	-18.4816	1.21159
+671.572	-18.5092	1.21241
+671.972	-18.5367	1.21323
+672.372	-18.5643	1.21405
+672.773	-18.5918	1.21485
+673.173	-18.6194	1.21565
+673.574	-18.6469	1.21644
+673.974	-18.6745	1.21722
+674.374	-18.702	1.218
+674.775	-18.7295	1.21877
+675.175	-18.7571	1.21954
+675.576	-18.7846	1.2203
+675.976	-18.8122	1.22105
+676.376	-18.8397	1.2218
+676.777	-18.8672	1.22254
+677.177	-18.8948	1.22328
+677.578	-18.9223	1.22401
+677.978	-18.9499	1.22474
+678.378	-18.9774	1.22546
+678.779	-19.0049	1.22618
+679.179	-19.0325	1.22689
+679.58	-19.06	1.2276
+679.98	-19.0875	1.22831
+680.38	-19.1151	1.22901
+680.781	-19.1426	1.22971
+681.181	-19.1701	1.2304
+681.582	-19.1977	1.2311
+681.982	-19.2252	1.23179
+682.382	-19.2527	1.23247
+682.783	-19.2803	1.23315
+683.183	-19.3078	1.23383
+683.584	-19.3353	1.23451
+683.984	-19.3629	1.23519
+684.384	-19.3904	1.23586
+684.785	-19.4179	1.23653
+685.185	-19.4455	1.2372
+685.586	-19.473	1.23787
+685.986	-19.5005	1.23854
+686.386	-19.528	1.2392
+686.787	-19.5556	1.23987
+687.187	-19.5831	1.24053
+687.588	-19.6106	1.24119
+687.988	-19.6382	1.24186
+688.388	-19.6657	1.24252
+688.789	-19.6932	1.24318

examples/Ag.txt

@@ -0,0 +1,244 @@
+#WL	Real	Imag
+0.1240	0.9980	0.0000
+0.1550	1.0000	0.0000
+0.2070	1.0040	0.0000
+0.2480	1.0020	0.0000
+0.2760	1.0020	0.0000
+0.3100	1.0060	0.0001
+0.3440	1.0060	0.0001
+0.3540	1.0060	0.0000
+0.4130	1.0020	0.0000
+0.4960	1.0020	0.0000
+0.6200	1.0020	0.0001
+0.8270	0.9980	0.0002
+1.2400	0.9960	0.0010
+1.5500	0.9940	0.0022
+1.6800	0.9960	0.0028
+2.0700	0.9940	0.0054
+2.4800	0.9940	0.0091
+3.1000	0.9959	0.0176
+3.3500	1.0140	0.0166
+3.5400	1.0020	0.0034
+4.1300	0.9920	0.0046
+4.4300	0.9880	0.0054
+4.7700	0.9841	0.0062
+5.1700	0.9781	0.0069
+5.6400	0.9742	0.0068
+6.2000	0.9643	0.0074
+6.8900	0.9565	0.0084
+7.2900	0.9467	0.0086
+7.7500	0.9370	0.0085
+8.2700	0.9274	0.0078
+8.8600	0.9120	0.0068
+9.5400	0.8892	0.0069
+10.3000	0.8574	0.0115
+11.3000	0.8133	0.0310
+11.8000	0.7914	0.0465
+12.4000	0.7659	0.0669
+13.1000	0.7381	0.0973
+13.8000	0.7140	0.1399
+14.6000	0.7034	0.1878
+15.5000	0.6998	0.2357
+15.9000	0.7020	0.2536
+16.3000	0.7026	0.2695
+16.8000	0.7001	0.2832
+17.2000	0.6914	0.3081
+17.7000	0.6991	0.3384
+18.2000	0.7154	0.3623
+18.8000	0.7352	0.3726
+19.4000	0.7391	0.3717
+20.0000	0.7316	0.3718
+20.7000	0.7133	0.3859
+21.4000	0.6973	0.4114
+22.1000	0.6916	0.4512
+23.0000	0.7081	0.4850
+23.4000	0.7165	0.4895
+23.8000	0.7141	0.4907
+24.3000	0.7060	0.4979
+24.8000	0.6974	0.5127
+25.3000	0.7432	0.3540
+25.8000	0.6931	0.5488
+26.4000	0.6970	0.5662
+27.0000	0.6985	0.5788
+28.2000	0.6966	0.6005
+29.5000	0.6828	0.6261
+31.0000	0.6674	0.6595
+31.8000	0.6581	0.6766
+32.6000	0.6451	0.6922
+33.5000	0.6232	0.7080
+34.4000	0.5927	0.7323
+35.4000	0.5421	0.7854
+36.5000	0.5335	0.8597
+37.6000	0.5440	0.9242
+38.7000	0.5695	0.9726
+40.0000	0.5835	0.9954
+41.3000	0.5741	1.0073
+42.8000	0.5343	1.0238
+43.5000	0.5027	1.0422
+44.3000	0.4655	1.0644
+45.1000	0.3447	1.0484
+45.9000	0.3625	1.1518
+45.9000	0.6187	0.7038
+46.8000	0.6078	0.7460
+47.7000	0.5896	0.7877
+48.6000	0.5749	0.8340
+50.6000	0.5484	0.9459
+52.8000	0.5482	1.1014
+53.9000	0.5834	1.2034
+55.1000	0.6695	1.2842
+56.4000	0.7774	1.3223
+57.7000	0.8778	1.2895
+59.0000	0.9152	1.2443
+60.5000	0.9196	1.2133
+62.0000	0.9053	1.2034
+65.3000	0.8820	1.2385
+68.9000	0.8606	1.3058
+72.9000	0.8534	1.4237
+77.5000	0.9049	1.5778
+80.0000	0.9673	1.6612
+82.7000	1.0648	1.7228
+85.5000	1.1799	1.7396
+88.6000	1.2805	1.7161
+91.8000	1.3533	1.6604
+95.4000	1.3841	1.5943
+99.2000	1.3835	1.5464
+102.5000	1.3733	1.5199
+105.1000	1.3663	1.5074
+107.8000	1.3617	1.4898
+110.7000	1.3515	1.4637
+113.7000	1.3248	1.4336
+117.0000	1.2866	1.4168
+120.4000	1.2494	1.4123
+124.0000	1.2175	1.4098
+127.8000	1.1901	1.3912
+134.8000	1.0946	1.3002
+137.8000	1.0155	1.2685
+140.9000	0.9196	1.2521
+144.2000	0.8138	1.2468
+147.6000	0.6929	1.2590
+151.2000	0.5596	1.2969
+155.0000	0.4270	1.3583
+159.0000	0.2907	1.4476
+163.1000	0.1820	1.5558
+167.5000	0.0804	1.6698
+172.2000	-0.0170	1.7917
+177.1000	-0.1119	1.9251
+182.3000	-0.2059	2.0737
+187.9000	-0.2869	2.2487
+193.7000	-0.3356	2.4261
+196.8000	-0.3658	2.5362
+200.0000	-0.3884	2.6586
+203.3000	-0.3820	2.7670
+206.6000	-0.3473	2.8575
+213.8000	-0.2882	3.0263
+221.4000	-0.2307	3.1408
+229.6000	-0.1835	3.2436
+238.4000	-0.1687	3.3649
+248.0000	-0.1377	3.5046
+253.0000	-0.0801	3.5640
+258.3000	-0.0189	3.6261
+263.8000	0.0599	3.7044
+269.5000	0.2023	3.7346
+275.5000	0.3604	3.7754
+281.8000	0.5910	3.7195
+288.3000	0.8399	3.5748
+295.2000	1.1410	3.2810
+298.8000	1.3324	3.0196
+302.4000	1.4601	2.6389
+306.1000	1.4639	2.1938
+310.0000	1.3317	1.7120
+311.5000	1.2091	1.4603
+313.9000	1.0286	1.2409
+315.5000	0.8257	1.0732
+317.9000	0.6146	0.9395
+319.5000	0.3875	0.8574
+322.0000	0.1820	0.8000
+323.7000	0.0086	0.7503
+326.3000	-0.1629	0.6975
+330.6000	-0.5233	0.6032
+332.4000	-0.7106	0.5791
+335.1000	-0.8858	0.5798
+339.7000	-1.1873	0.5802
+344.4000	-1.4810	0.5902
+354.2000	-2.0299	0.6019
+364.7000	-2.5575	0.5989
+375.7000	-2.7489	0.6680
+387.5000	-3.2392	0.6950
+400.0000	-3.7726	0.6747
+413.3000	-4.4222	0.7301
+427.5000	-5.0820	0.7232
+442.8000	-5.7354	0.7536
+459.2000	-6.5329	0.7373
+476.9000	-7.3810	0.7181
+495.9000	-8.2775	0.7488
+516.6000	-9.4080	0.7982
+539.1000	-10.5459	0.8385
+563.6000	-11.8881	0.8280
+590.4000	-13.3810	0.8857
+619.9000	-15.0372	1.0166
+652.6000	-17.2029	1.1620
+688.8000	-19.6940	1.2432
+729.3000	-22.4457	1.4030
+774.9000	-25.8877	1.4557
+826.6000	-30.2290	1.5950
+885.6000	-35.3759	1.9397
+953.7000	-41.3057	2.5463
+1033.0000	-48.8090	3.1595
+1127.0000	-58.7659	3.8503
+1240.0000	-71.9719	5.5864
+1265.0000	-60.3878	5.8350
+1291.0000	-62.5797	6.0667
+1305.0000	-79.9743	6.4082
+1319.0000	-64.8099	6.3190
+1348.0000	-67.2433	6.5844
+1378.0000	-69.8880	6.8801
+1378.0000	-89.7136	7.5082
+1409.0000	-69.8797	7.0475
+1442.0000	-75.5042	7.4994
+1459.0000	-101.8111	9.0092
+1476.0000	-78.6590	7.8499
+1512.0000	-82.2394	8.2628
+1550.0000	-86.6424	8.7422
+1550.0000	-116.3758	11.1024
+1590.0000	-91.3497	9.2829
+1631.0000	-96.5746	9.8597
+1653.0000	-131.8606	14.3520
+1675.0000	-101.7406	10.4838
+1722.0000	-107.8716	11.1696
+1771.0000	-114.1798	11.9198
+1771.0000	-148.1277	20.5936
+1823.0000	-122.8759	12.8316
+1879.0000	-129.6000	13.6800
+1907.0000	-176.1279	23.2218
+1937.0000	-138.8506	14.7264
+2000.0000	-148.4175	15.8600
+2066.0000	-158.3138	16.8336
+2066.0000	-206.2279	30.6432
+2138.0000	-168.4686	18.9540
+2214.0000	-181.6509	20.8980
+2296.0000	-195.3227	23.0440
+2384.0000	-209.4791	25.4620
+2480.0000	-227.1283	28.3578
+2583.0000	-245.4760	31.6198
+2695.0000	-267.7871	35.5224
+2818.0000	-291.0458	39.9456
+2952.0000	-318.8098	45.2870
+3100.0000	-351.5162	52.1512
+3263.0000	-389.6807	60.8256
+3444.0000	-433.8859	71.4780
+3647.0000	-484.7428	84.6430
+3875.0000	-547.5844	101.5200
+4133.0000	-624.0271	122.7892
+4428.0000	-715.8482	149.8868
+4769.0000	-830.7472	185.7160
+5166.0000	-965.7622	233.6232
+5636.0000	-1136.4194	300.9000
+6199.0000	-1340.3240	396.2700
+6526.0000	-1454.4384	460.1120
+6888.0000	-1587.6711	538.9360
+7293.0000	-1750.5835	634.1850
+7749.0000	-1936.8826	750.4896
+8266.0000	-2129.2775	889.3422
+8856.0000	-2326.0839	1056.1720
+9537.0000	-2575.7559	1274.7240
+9919.0000	-2711.8179	1408.0140

examples/Si-int.txt

@@ -0,0 +1,988 @@
+301.201	8.62319	40.7522
+301.602	8.88308	40.48
+302.002	9.13357	40.2111
+302.402	9.36932	39.9476
+302.803	9.58666	39.6908
+303.203	9.78861	39.4398
+303.604	9.97983	39.1934
+304.004	10.165	38.9502
+304.404	10.3467	38.7097
+304.805	10.5238	38.4732
+305.205	10.695	38.242
+305.606	10.8586	38.0176
+306.006	11.0147	37.8
+306.406	11.1644	37.5883
+306.807	11.3089	37.3815
+307.207	11.4493	37.1787
+307.608	11.5844	36.9807
+308.008	11.7121	36.7892
+308.408	11.8301	36.6057
+308.809	11.9375	36.4313
+309.209	12.0379	36.2636
+309.61	12.1359	36.1
+310.01	12.2361	35.9377
+310.41	12.3418	35.7749
+310.811	12.4508	35.6143
+311.211	12.5597	35.4595
+311.612	12.6651	35.3141
+312.012	12.7648	35.1793
+312.412	12.8592	35.0518
+312.813	12.9489	34.928
+313.213	13.0343	34.8044
+313.614	13.1164	34.6794
+314.014	13.1963	34.5553
+314.414	13.2754	34.4344
+314.815	13.355	34.3195
+315.215	13.4357	34.2115
+315.616	13.5169	34.109
+316.016	13.598	34.0101
+316.416	13.6784	33.9132
+316.817	13.7575	33.8171
+317.217	13.8353	33.7226
+317.618	13.9116	33.6304
+318.018	13.9864	33.5411
+318.418	14.0598	33.4552
+318.819	14.1321	33.3726
+319.219	14.204	33.293
+319.62	14.2758	33.2164
+320.02	14.348	33.1426
+320.42	14.4208	33.0715
+320.821	14.4945	33.0028
+321.221	14.5692	32.9365
+321.622	14.6452	32.8724
+322.022	14.722	32.8102
+322.422	14.7992	32.7497
+322.823	14.8762	32.6906
+323.223	14.9528	32.6327
+323.624	15.0288	32.5761
+324.024	15.1042	32.5208
+324.424	15.1793	32.467
+324.825	15.2542	32.4149
+325.225	15.3288	32.3643
+325.626	15.4033	32.3149
+326.026	15.478	32.2664
+326.426	15.5527	32.2186
+326.827	15.6277	32.1713
+327.227	15.703	32.1248
+327.628	15.7783	32.0793
+328.028	15.8538	32.0353
+328.428	15.9295	31.993
+328.829	16.0056	31.9528
+329.229	16.0826	31.915
+329.63	16.1608	31.8801
+330.03	16.2404	31.8483
+330.43	16.3209	31.8189
+330.831	16.4016	31.7908
+331.231	16.4817	31.7628
+331.632	16.5602	31.7341
+332.032	16.6371	31.7042
+332.432	16.713	31.6743
+332.833	16.7886	31.6456
+333.233	16.8649	31.6193
+333.634	16.9424	31.5962
+334.034	17.021	31.5764
+334.434	17.1004	31.5591
+334.835	17.1802	31.5436
+335.235	17.2598	31.5294
+335.636	17.3392	31.5159
+336.036	17.4181	31.5032
+336.436	17.4965	31.4911
+336.837	17.5743	31.4797
+337.237	17.6516	31.469
+337.638	17.7286	31.4593
+338.038	17.8059	31.4513
+338.438	17.8838	31.4454
+338.839	17.9629	31.4424
+339.239	18.0434	31.4425
+339.64	18.1249	31.4452
+340.04	18.2069	31.45
+340.44	18.2888	31.4564
+340.841	18.3702	31.4638
+341.241	18.451	31.4723
+341.642	18.5319	31.4825
+342.042	18.6135	31.4952
+342.442	18.6961	31.5109
+342.843	18.7803	31.5301
+343.243	18.8664	31.5529
+343.644	18.9543	31.5793
+344.044	19.044	31.6091
+344.444	19.1357	31.6421
+344.845	19.2293	31.6783
+345.245	19.324	31.7175
+345.646	19.4194	31.7597
+346.046	19.5146	31.8046
+346.446	19.609	31.8523
+346.847	19.7031	31.903
+347.247	19.7987	31.9575
+347.648	19.8979	32.0168
+348.048	20.0027	32.0817
+348.448	20.1151	32.153
+348.849	20.2356	32.231
+349.249	20.3628	32.3148
+349.65	20.4948	32.4037
+350.05	20.6301	32.4966
+350.45	20.7671	32.5929
+350.851	20.9064	32.693
+351.251	21.0501	32.7985
+351.652	21.2005	32.911
+352.052	21.3597	33.0319
+352.452	21.5298	33.1626
+352.853	21.711	33.3026
+353.253	21.9023	33.4501
+353.654	22.1024	33.603
+354.054	22.3104	33.7595
+354.454	22.5253	33.9178
+354.855	22.7496	34.078
+355.255	22.988	34.2417
+355.656	23.2452	34.4104
+356.056	23.5259	34.5855
+356.456	23.835	34.7687
+356.857	24.1739	34.9589
+357.257	24.54	35.1519
+357.658	24.9305	35.3431
+358.058	25.3425	35.5282
+358.458	25.7732	35.7026
+358.859	26.2224	35.863
+359.259	26.6932	36.0075
+359.66	27.1891	36.1341
+360.06	27.7135	36.2409
+360.46	28.2698	36.3262
+360.861	28.8589	36.3879
+361.261	29.4756	36.4238
+361.662	30.1137	36.4316
+362.062	30.7667	36.409
+362.462	31.4284	36.3536
+362.863	32.0936	36.2631
+363.263	32.7624	36.1335
+363.664	33.4365	35.9608
+364.064	34.1178	35.7408
+364.464	34.8083	35.4694
+364.865	35.5096	35.1425
+365.265	36.2176	34.7602
+365.666	36.9214	34.3283
+366.066	37.6094	33.8524
+366.466	38.27	33.3386
+366.867	38.8917	32.7927
+367.267	39.4668	32.2195
+367.668	39.9972	31.6217
+368.068	40.486	31.0012
+368.468	40.9365	30.3601
+368.869	41.352	29.7007
+369.269	41.7356	29.0248
+369.67	42.088	28.3331
+370.07	42.4087	27.6255
+370.47	42.6972	26.9017
+370.871	42.953	26.1614
+371.271	43.1755	25.4045
+371.672	43.3639	24.6334
+372.072	43.5172	23.8568
+372.472	43.634	23.0841
+372.873	43.7133	22.325
+373.273	43.7538	21.589
+373.674	43.7546	20.885
+374.074	43.7178	20.2134
+374.474	43.6474	19.5698
+374.875	43.5475	18.9495
+375.275	43.4221	18.348
+375.676	43.2752	17.7605
+376.076	43.1105	17.1836
+376.476	42.9301	16.6182
+376.877	42.7358	16.0671
+377.277	42.5291	15.5327
+377.678	42.3117	15.0176
+378.078	42.0854	14.5244
+378.478	41.8516	14.0546
+378.879	41.6118	13.6078
+379.279	41.3677	13.183
+379.68	41.1206	12.7795
+380.08	40.8719	12.3964
+380.48	40.6233	12.0327
+380.881	40.3749	11.6873
+381.281	40.1261	11.3581
+381.682	39.8759	11.0434
+382.082	39.6235	10.7412
+382.482	39.368	10.4497
+382.883	39.1086	10.1669
+383.283	38.8461	9.89242
+383.684	38.5828	9.62711
+384.084	38.3214	9.372
+384.484	38.0645	9.12809
+384.885	37.8146	8.89639
+385.285	37.5741	8.67774
+385.686	37.3424	8.47127
+386.086	37.1178	8.27515
+386.486	36.8981	8.08753
+386.887	36.6814	7.90655
+387.287	36.4658	7.73039
+387.688	36.2492	7.55723
+388.088	36.0313	7.38661
+388.488	35.8134	7.21949
+388.889	35.5967	7.05692
+389.289	35.3826	6.89995
+389.69	35.1724	6.74961
+390.09	34.9674	6.60691
+390.49	34.768	6.47188
+390.891	34.5735	6.34347
+391.291	34.3834	6.22061
+391.692	34.1969	6.10223
+392.092	34.0134	5.98723
+392.492	33.8323	5.87454
+392.893	33.653	5.76349
+393.293	33.4757	5.65421
+393.694	33.3004	5.54696
+394.094	33.1272	5.44197
+394.494	32.9563	5.3395
+394.895	32.7877	5.23978
+395.295	32.6216	5.14311
+395.696	32.4581	5.04993
+396.096	32.2973	4.96075
+396.496	32.1395	4.87607
+396.897	31.9846	4.79639
+397.297	31.833	4.72221
+397.698	31.6847	4.65389
+398.098	31.5394	4.59056
+398.498	31.3966	4.5307
+398.899	31.2559	4.4728
+399.299	31.1166	4.41532
+399.7	30.9784	4.35675
+400.1	30.8407	4.29556
+400.501	30.7034	4.23111
+400.901	30.5672	4.16443
+401.301	30.4326	4.09672
+401.702	30.3005	4.02921
+402.102	30.1715	3.96311
+402.503	30.0463	3.89964
+402.903	29.9253	3.83981
+403.303	29.8078	3.78372
+403.704	29.6928	3.7313
+404.104	29.5793	3.68246
+404.505	29.4663	3.63711
+404.905	29.3527	3.59515
+405.305	29.2377	3.55651
+405.706	29.1208	3.52083
+406.106	29.0028	3.48728
+406.507	28.8847	3.45496
+406.907	28.7675	3.423
+407.307	28.652	3.3905
+407.708	28.5393	3.35659
+408.108	28.4302	3.32046
+408.509	28.3246	3.28218
+408.909	28.222	3.24218
+409.309	28.1221	3.20093
+409.71	28.0244	3.15888
+410.11	27.9283	3.11648
+410.511	27.8335	3.07418
+410.911	27.7395	3.03242
+411.311	27.6464	2.99147
+411.712	27.5541	2.95162
+412.112	27.4626	2.91311
+412.513	27.3721	2.87623
+412.913	27.2826	2.84123
+413.313	27.194	2.80838
+413.714	27.1064	2.77786
+414.114	27.0198	2.74952
+414.515	26.934	2.72314
+414.915	26.8491	2.6985
+415.315	26.7649	2.67538
+415.716	26.6814	2.65356
+416.116	26.5985	2.63282
+416.517	26.5163	2.61293
+416.917	26.4347	2.59354
+417.317	26.3537	2.57432
+417.718	26.2736	2.55491
+418.118	26.1944	2.53497
+418.519	26.116	2.51416
+418.919	26.0387	2.49212
+419.319	25.9625	2.46861
+419.72	25.8873	2.44373
+420.12	25.8132	2.41765
+420.521	25.7402	2.39055
+420.921	25.6682	2.3626
+421.321	25.5974	2.33398
+421.722	25.5276	2.30486
+422.122	25.4589	2.27547
+422.523	25.3911	2.24622
+422.923	25.3239	2.21755
+423.323	25.2573	2.18991
+423.724	25.1909	2.16375
+424.124	25.1245	2.1395
+424.525	25.0581	2.11762
+424.925	24.9914	2.09847
+425.325	24.9245	2.08182
+425.726	24.8578	2.06719
+426.126	24.7915	2.05409
+426.527	24.7258	2.04205
+426.927	24.661	2.03057
+427.327	24.5973	2.01918
+427.728	24.5351	2.00739
+428.128	24.4743	1.995
+428.529	24.4146	1.98195
+428.929	24.3559	1.96823
+429.329	24.2978	1.95381
+429.73	24.2401	1.93867
+430.13	24.1826	1.92277
+430.531	24.125	1.90611
+430.931	24.0672	1.88878
+431.331	24.0093	1.87135
+431.732	23.9514	1.85446
+432.132	23.8939	1.83874
+432.533	23.8368	1.82483
+432.933	23.7805	1.81339
+433.333	23.725	1.80505
+433.734	23.6705	1.80038
+434.134	23.6172	1.79897
+434.535	23.5647	1.79959
+434.935	23.5131	1.80098
+435.335	23.4622	1.80189
+435.736	23.4118	1.80107
+436.136	23.3619	1.79726
+436.537	23.3122	1.78922
+436.937	23.2627	1.77594
+437.337	23.2133	1.75807
+437.738	23.1641	1.73691
+438.138	23.115	1.71376
+438.539	23.066	1.68993
+438.939	23.0172	1.66671
+439.339	22.9685	1.64541
+439.74	22.92	1.62733
+440.14	22.8716	1.61324
+440.541	22.8235	1.60235
+440.941	22.7758	1.59359
+441.341	22.7287	1.58588
+441.742	22.6824	1.57814
+442.142	22.6369	1.56928
+442.543	22.5924	1.55824
+442.943	22.5491	1.54394
+443.343	22.5069	1.52621
+443.744	22.4657	1.50613
+444.144	22.425	1.48491
+444.545	22.3847	1.46372
+444.945	22.3446	1.44378
+445.345	22.3043	1.42626
+445.746	22.2635	1.41236
+446.146	22.2221	1.40326
+446.547	22.18	1.3991
+446.947	22.1373	1.39856
+447.347	22.0944	1.40021
+447.748	22.0514	1.40262
+448.148	22.0086	1.40436
+448.549	21.9663	1.404
+448.949	21.9246	1.40011
+449.349	21.8838	1.39129
+449.75	21.8441	1.37719
+450.15	21.8052	1.35897
+450.551	21.767	1.33792
+450.951	21.7294	1.31529
+451.351	21.6922	1.29236
+451.752	21.6552	1.2704
+452.152	21.6183	1.25069
+452.553	21.5814	1.23449
+452.953	21.5442	1.2226
+453.353	21.507	1.21456
+453.754	21.4697	1.20968
+454.154	21.4325	1.2073
+454.555	21.3953	1.20673
+454.955	21.3583	1.20729
+455.355	21.3216	1.2083
+455.756	21.2851	1.20909
+456.156	21.2489	1.20906
+456.557	21.2132	1.20819
+456.957	21.1778	1.20661
+457.357	21.1429	1.20447
+457.758	21.1083	1.20192
+458.158	21.0741	1.1991
+458.559	21.0403	1.19617
+458.959	21.0069	1.19326
+459.359	20.9739	1.19052
+459.76	20.9413	1.18811
+460.16	20.909	1.18616
+460.561	20.8769	1.18482
+460.961	20.8448	1.18424
+461.361	20.8128	1.18457
+461.762	20.7807	1.18595
+462.162	20.7483	1.18853
+462.563	20.7157	1.19245
+462.963	20.6828	1.19775
+463.363	20.6496	1.2039
+463.764	20.6163	1.21019
+464.164	20.5832	1.21588
+464.565	20.5504	1.22026
+464.965	20.5181	1.22261
+465.365	20.4864	1.2222
+465.766	20.4555	1.21832
+466.166	20.4257	1.21024
+466.567	20.3969	1.19776
+466.967	20.369	1.18188
+467.367	20.3418	1.16381
+467.768	20.3149	1.14477
+468.168	20.2882	1.12594
+468.569	20.2615	1.10855
+468.969	20.2344	1.09379
+469.369	20.2069	1.08288
+469.77	20.1785	1.07697
+470.17	20.1493	1.07605
+470.571	20.1197	1.07865
+470.971	20.0897	1.08319
+471.371	20.0598	1.0881
+471.772	20.0301	1.09184
+472.172	20.0011	1.09282
+472.573	19.9728	1.08949
+472.973	19.9456	1.08027
+473.373	19.9198	1.06366
+473.774	19.8953	1.0396
+474.174	19.872	1.00986
+474.575	19.8496	0.97625
+474.975	19.8278	0.940629
+475.375	19.8063	0.904834
+475.776	19.7849	0.870705
+476.176	19.7634	0.840084
+476.577	19.7415	0.814811
+476.977	19.7188	0.796724
+477.377	19.6954	0.786734
+477.778	19.6712	0.783646
+478.178	19.6464	0.785969
+478.579	19.6212	0.792214
+478.979	19.5957	0.800893
+479.379	19.57	0.810515
+479.78	19.5443	0.819592
+480.18	19.5187	0.826634
+480.581	19.4934	0.830153
+480.981	19.4685	0.828988
+481.381	19.444	0.823508
+481.782	19.4198	0.814522
+482.182	19.396	0.80284
+482.583	19.3724	0.78927
+482.983	19.3491	0.774623
+483.383	19.3261	0.759707
+483.784	19.3032	0.745332
+484.184	19.2805	0.732307
+484.585	19.258	0.721377
+484.985	19.2356	0.712645
+485.385	19.2133	0.705847
+485.786	19.1911	0.700713
+486.186	19.169	0.696975
+486.587	19.147	0.694363
+486.987	19.1251	0.692609
+487.387	19.1032	0.691443
+487.788	19.0814	0.690597
+488.188	19.0597	0.689802
+488.589	19.0379	0.688858
+488.989	19.0163	0.68771
+489.389	18.9947	0.686322
+489.79	18.9731	0.684656
+490.19	18.9516	0.682677
+490.591	18.9302	0.680347
+490.991	18.9089	0.67763
+491.391	18.8877	0.674489
+491.792	18.8666	0.670889
+492.192	18.8457	0.666794
+492.593	18.8248	0.662271
+492.993	18.8041	0.657486
+493.393	18.7835	0.652609
+493.794	18.763	0.647812
+494.194	18.7426	0.643266
+494.595	18.7223	0.639143
+494.995	18.7022	0.635614
+495.395	18.6821	0.632849
+495.796	18.6622	0.631021
+496.196	18.6424	0.630268
+496.597	18.6228	0.630441
+496.997	18.6033	0.631236
+497.397	18.5839	0.632348
+497.798	18.5647	0.633473
+498.198	18.5457	0.634306
+498.599	18.5269	0.634543
+498.999	18.5084	0.633878
+499.399	18.49	0.632007
+499.8	18.472	0.628626
+500.2	18.4541	0.623484
+500.601	18.4366	0.616806
+501.001	18.4191	0.609062
+501.401	18.4018	0.600723
+501.802	18.3846	0.592259
+502.202	18.3673	0.584144
+502.603	18.35	0.576848
+503.003	18.3326	0.570843
+503.403	18.315	0.5666
+503.804	18.2972	0.56459
+504.204	18.2791	0.565262
+504.605	18.2607	0.568527
+505.005	18.2421	0.573781
+505.405	18.2234	0.580394
+505.806	18.2047	0.587738
+506.206	18.1859	0.595186
+506.607	18.1673	0.602109
+507.007	18.1488	0.607879
+507.407	18.1306	0.611868
+507.808	18.1127	0.613448
+508.208	18.0951	0.611995
+508.609	18.0781	0.607281
+509.009	18.0613	0.599804
+509.409	18.045	0.590142
+509.81	18.0289	0.578873
+510.21	18.013	0.566573
+510.611	17.9973	0.55382
+511.011	17.9817	0.541191
+511.411	17.9661	0.529264
+511.812	17.9506	0.518615
+512.212	17.935	0.509823
+512.613	17.9194	0.503389
+513.013	17.9037	0.499236
+513.413	17.8879	0.497008
+513.814	17.872	0.496351
+514.214	17.8562	0.496907
+514.615	17.8403	0.49832
+515.015	17.8245	0.500235
+515.415	17.8087	0.502296
+515.816	17.793	0.504147
+516.216	17.7774	0.505431
+516.617	17.762	0.505792
+517.017	17.7467	0.504977
+517.417	17.7315	0.503095
+517.818	17.7164	0.500332
+518.218	17.7015	0.496877
+518.619	17.6866	0.492919
+519.019	17.6718	0.488644
+519.419	17.6572	0.484241
+519.82	17.6425	0.479899
+520.22	17.6279	0.475805
+520.621	17.6134	0.472146
+521.021	17.5989	0.46911
+521.421	17.5844	0.466755
+521.822	17.57	0.464921
+522.222	17.5555	0.463431
+522.623	17.5412	0.462106
+523.023	17.5269	0.460765
+523.423	17.5128	0.459232
+523.824	17.4987	0.457325
+524.224	17.4847	0.454868
+524.625	17.4709	0.45168
+525.025	17.4572	0.447582
+525.425	17.4437	0.442397
+525.826	17.4303	0.43603
+526.226	17.417	0.428669
+526.627	17.4039	0.420562
+527.027	17.3908	0.411955
+527.427	17.3777	0.403094
+527.828	17.3647	0.394228
+528.228	17.3515	0.385602
+528.629	17.3383	0.377464
+529.029	17.325	0.370059
+529.429	17.3115	0.363636
+529.83	17.2979	0.35844
+530.23	17.284	0.354682
+530.631	17.2699	0.352325
+531.031	17.2557	0.35123
+531.431	17.2413	0.351254
+531.832	17.227	0.352258
+532.232	17.2127	0.354101
+532.633	17.1984	0.356641
+533.033	17.1843	0.359739
+533.433	17.1703	0.363253
+533.834	17.1565	0.367043
+534.234	17.1431	0.370969
+534.635	17.1299	0.37489
+535.035	17.1171	0.378709
+535.435	17.1045	0.382355
+535.836	17.0922	0.385758
+536.236	17.0801	0.388852
+536.637	17.0682	0.391566
+537.037	17.0565	0.393832
+537.437	17.0448	0.395581
+537.838	17.0333	0.396745
+538.238	17.0217	0.397254
+538.639	17.0102	0.39704
+539.039	16.9986	0.396035
+539.439	16.9869	0.394192
+539.84	16.9752	0.391617
+540.24	16.9635	0.38847
+540.641	16.9516	0.384915
+541.041	16.9397	0.381114
+541.441	16.9278	0.37723
+541.842	16.9159	0.373424
+542.242	16.9039	0.369861
+542.643	16.8919	0.366701
+543.043	16.8799	0.364107
+543.443	16.868	0.362243
+543.844	16.856	0.36127
+544.244	16.844	0.361259
+544.645	16.8321	0.362018
+545.045	16.8202	0.363307
+545.445	16.8084	0.364886
+545.846	16.7966	0.366515
+546.246	16.7849	0.367954
+546.647	16.7733	0.368963
+547.047	16.7617	0.369302
+547.447	16.7503	0.368731
+547.848	16.7389	0.367011
+548.248	16.7277	0.363901
+548.649	16.7166	0.359161
+549.049	16.7057	0.352689
+549.449	16.6949	0.344761
+549.85	16.6841	0.335721
+550.25	16.6735	0.325912
+550.651	16.6629	0.315677
+551.051	16.6524	0.305359
+551.451	16.6419	0.295303
+551.852	16.6315	0.285851
+552.252	16.621	0.277346
+552.653	16.6105	0.270131
+553.053	16.6001	0.26455
+553.453	16.5895	0.260947
+553.854	16.5789	0.259586
+554.254	16.5683	0.260289
+554.655	16.5577	0.262723
+555.055	16.547	0.266553
+555.455	16.5364	0.271444
+555.856	16.5258	0.277062
+556.256	16.5152	0.283072
+556.657	16.5047	0.289141
+557.057	16.4943	0.294933
+557.457	16.484	0.300114
+557.858	16.4738	0.30435
+558.258	16.4638	0.307307
+558.659	16.4539	0.308655
+559.059	16.4442	0.308304
+559.459	16.4346	0.306472
+559.86	16.4252	0.3034
+560.26	16.4158	0.299329
+560.661	16.4065	0.294497
+561.061	16.3971	0.289147
+561.461	16.3878	0.283518
+561.862	16.3785	0.27785
+562.262	16.369	0.272384
+562.663	16.3595	0.26736
+563.063	16.3499	0.263019
+563.463	16.3401	0.2596
+563.864	16.3302	0.257325
+564.264	16.32	0.256189
+564.665	16.3098	0.256035
+565.065	16.2994	0.256705
+565.465	16.2891	0.258042
+565.866	16.2786	0.259886
+566.266	16.2682	0.262079
+566.667	16.2579	0.264465
+567.067	16.2476	0.266883
+567.467	16.2375	0.269176
+567.868	16.2276	0.271186
+568.268	16.2178	0.272754
+568.669	16.2083	0.273722
+569.069	16.199	0.273962
+569.469	16.19	0.273496
+569.87	16.1812	0.272395
+570.27	16.1725	0.270728
+570.671	16.164	0.268565
+571.071	16.1557	0.265977
+571.471	16.1474	0.263033
+571.872	16.1392	0.259804
+572.272	16.131	0.256359
+572.673	16.1228	0.25277
+573.073	16.1146	0.249105
+573.473	16.1063	0.245435
+573.874	16.0979	0.241831
+574.274	16.0894	0.238359
+574.675	16.0807	0.235057
+575.075	16.072	0.231945
+575.475	16.0631	0.229044
+575.876	16.0542	0.226373
+576.276	16.0452	0.223951
+576.677	16.0362	0.221799
+577.077	16.0272	0.219935
+577.477	16.0181	0.218381
+577.878	16.0091	0.217154
+578.278	16.0001	0.216276
+578.679	15.9912	0.215766
+579.079	15.9823	0.215643
+579.479	15.9735	0.215928
+579.88	15.9649	0.216619
+580.28	15.9563	0.217673
+580.681	15.9478	0.219037
+581.081	15.9394	0.220661
+581.481	15.9311	0.222494
+581.882	15.9228	0.224484
+582.282	15.9145	0.226581
+582.683	15.9063	0.228734
+583.083	15.8982	0.230891
+583.483	15.8901	0.233002
+583.884	15.8819	0.235015
+584.284	15.8738	0.23688
+584.685	15.8657	0.238545
+585.085	15.8576	0.239963
+585.485	15.8495	0.241122
+585.886	15.8414	0.242034
+586.286	15.8333	0.242706
+586.687	15.8252	0.24315
+587.087	15.8171	0.243375
+587.487	15.809	0.243391
+587.888	15.801	0.243208
+588.288	15.7931	0.242836
+588.689	15.7851	0.242285
+589.089	15.7773	0.241564
+589.489	15.7695	0.240683
+589.89	15.7618	0.239653
+590.29	15.7541	0.238484
+590.691	15.7466	0.237184
+591.091	15.7391	0.235768
+591.491	15.7317	0.234247
+591.892	15.7244	0.232636
+592.292	15.7171	0.230948
+592.693	15.7099	0.229196
+593.093	15.7027	0.227394
+593.493	15.6956	0.225556
+593.894	15.6884	0.223694
+594.294	15.6813	0.221823
+594.695	15.6742	0.219955
+595.095	15.6671	0.218104
+595.495	15.66	0.216284
+595.896	15.6529	0.214508
+596.296	15.6457	0.212789
+596.697	15.6385	0.211136
+597.097	15.6313	0.209553
+597.497	15.6241	0.208041
+597.898	15.6168	0.206606
+598.298	15.6096	0.20525
+598.699	15.6024	0.203976
+599.099	15.5953	0.202788
+599.499	15.5881	0.201689
+599.9	15.581	0.200682
+600.3	15.574	0.199771
+600.701	15.567	0.198959
+601.101	15.5601	0.198248
+601.502	15.5533	0.197644
+601.902	15.5466	0.197148
+602.302	15.54	0.196762
+602.703	15.5334	0.196475
+603.103	15.527	0.196278
+603.504	15.5205	0.196157
+603.904	15.5142	0.196101
+604.304	15.5079	0.196098
+604.705	15.5016	0.196137
+605.105	15.4953	0.196205
+605.506	15.489	0.196292
+605.906	15.4827	0.196385
+606.306	15.4763	0.196472
+606.707	15.4699	0.196542
+607.107	15.4635	0.196582
+607.508	15.4569	0.196582
+607.908	15.4503	0.19653
+608.308	15.4437	0.196417
+608.709	15.4369	0.196241
+609.109	15.4301	0.196002
+609.51	15.4232	0.195697
+609.91	15.4163	0.195327
+610.31	15.4094	0.19489
+610.711	15.4024	0.194384
+611.111	15.3955	0.193809
+611.512	15.3886	0.193163
+611.912	15.3817	0.192446
+612.312	15.3749	0.191656
+612.713	15.3681	0.190792
+613.113	15.3614	0.189853
+613.514	15.3548	0.188838
+613.914	15.3483	0.187745
+614.314	15.3419	0.18658
+614.715	15.3356	0.185358
+615.115	15.3294	0.184096
+615.516	15.3232	0.182811
+615.916	15.3171	0.181518
+616.316	15.311	0.180234
+616.717	15.305	0.178977
+617.117	15.299	0.177762
+617.518	15.2929	0.176606
+617.918	15.2869	0.175527
+618.318	15.2808	0.174539
+618.719	15.2747	0.173661
+619.119	15.2686	0.172908
+619.52	15.2623	0.172297
+619.92	15.256	0.171846
+620.32	15.2496	0.171561
+620.721	15.2432	0.171424
+621.121	15.2366	0.171411
+621.522	15.2301	0.171496
+621.922	15.2235	0.171654
+622.322	15.2168	0.17186
+622.723	15.2102	0.172089
+623.123	15.2036	0.172316
+623.524	15.197	0.172515
+623.924	15.1905	0.172663
+624.324	15.184	0.172733
+624.725	15.1776	0.172701
+625.125	15.1713	0.172542
+625.526	15.1651	0.172231
+625.926	15.159	0.171742
+626.326	15.1531	0.171051
+626.727	15.1473	0.17015
+627.127	15.1416	0.169061
+627.528	15.1361	0.167806
+627.928	15.1306	0.16641
+628.328	15.1252	0.164895
+628.729	15.12	0.163284
+629.129	15.1147	0.161603
+629.53	15.1095	0.159873
+629.93	15.1043	0.158118
+630.33	15.0992	0.156362
+630.731	15.094	0.154629
+631.131	15.0889	0.15294
+631.532	15.0837	0.151321
+631.932	15.0784	0.149794
+632.332	15.0731	0.148383
+632.733	15.0678	0.147111
+633.133	15.0623	0.145988
+633.534	15.0568	0.145005
+633.934	15.0513	0.144149
+634.334	15.0456	0.143409
+634.735	15.0399	0.142773
+635.135	15.0342	0.142229
+635.536	15.0284	0.141766
+635.936	15.0227	0.141372
+636.336	15.0168	0.141035
+636.737	15.011	0.140744
+637.137	15.0052	0.140486
+637.538	14.9993	0.14025
+637.938	14.9935	0.140025
+638.338	14.9876	0.139798
+638.739	14.9818	0.139558
+639.139	14.976	0.139293
+639.54	14.9702	0.138994
+639.94	14.9645	0.138662
+640.34	14.9588	0.138298
+640.741	14.9531	0.137905
+641.141	14.9475	0.137484
+641.542	14.9418	0.137037
+641.942	14.9363	0.136566
+642.342	14.9307	0.136073
+642.743	14.9252	0.135559
+643.143	14.9197	0.135028
+643.544	14.9142	0.134479
+643.944	14.9087	0.133917
+644.344	14.9033	0.133341
+644.745	14.8979	0.132755
+645.145	14.8926	0.13216
+645.546	14.8872	0.131558
+645.946	14.8819	0.13095
+646.346	14.8767	0.13034
+646.747	14.8714	0.129731
+647.147	14.8662	0.129126
+647.548	14.8611	0.128529
+647.948	14.8559	0.127942
+648.348	14.8508	0.127371
+648.749	14.8458	0.126818
+649.149	14.8407	0.126286
+649.55	14.8357	0.125779
+649.95	14.8308	0.1253
+650.35	14.8259	0.124854
+650.751	14.821	0.124443
+651.151	14.8162	0.124071
+651.552	14.8114	0.123741
+651.952	14.8067	0.123458
+652.352	14.802	0.123223
+652.753	14.7974	0.123041
+653.153	14.7928	0.122911
+653.554	14.7883	0.122827
+653.954	14.7838	0.122783
+654.354	14.7793	0.122772
+654.755	14.7748	0.122788
+655.155	14.7704	0.122825
+655.556	14.766	0.122876
+655.956	14.7616	0.122934
+656.356	14.7572	0.122994
+656.757	14.7528	0.123049
+657.157	14.7484	0.123093
+657.558	14.744	0.123119
+657.958	14.7396	0.12312
+658.358	14.7351	0.123091
+658.759	14.7307	0.123025
+659.159	14.7262	0.122916
+659.56	14.7216	0.122757
+659.96	14.7171	0.122545
+660.36	14.7125	0.122282
+660.761	14.7078	0.121973
+661.161	14.7032	0.121621
+661.562	14.6985	0.121232
+661.962	14.6937	0.120808
+662.362	14.689	0.120355
+662.763	14.6842	0.119876
+663.163	14.6794	0.119375
+663.564	14.6747	0.118857
+663.964	14.6698	0.118325
+664.364	14.665	0.117783
+664.765	14.6602	0.117237
+665.165	14.6554	0.116689
+665.566	14.6505	0.116145
+665.966	14.6457	0.115607
+666.366	14.6409	0.11508
+666.767	14.636	0.114569
+667.167	14.6312	0.114075
+667.568	14.6264	0.113595
+667.968	14.6216	0.11313
+668.368	14.6168	0.112676
+668.769	14.6121	0.112232
+669.169	14.6073	0.111797
+669.57	14.6026	0.111369
+669.97	14.5979	0.110945
+670.37	14.5932	0.110525
+670.771	14.5885	0.110107
+671.171	14.5839	0.109688
+671.572	14.5792	0.109268
+671.972	14.5746	0.108845
+672.372	14.5701	0.108416
+672.773	14.5655	0.10798
+673.173	14.561	0.107536
+673.574	14.5565	0.107082
+673.974	14.5521	0.106616
+674.374	14.5477	0.106138
+674.775	14.5433	0.105652
+675.175	14.539	0.105159
+675.576	14.5347	0.104663
+675.976	14.5304	0.104166
+676.376	14.5262	0.103672
+676.777	14.522	0.103183
+677.177	14.5178	0.102702
+677.578	14.5137	0.102231
+677.978	14.5096	0.101774
+678.378	14.5055	0.101334
+678.779	14.5015	0.100912
+679.179	14.4975	0.100513
+679.58	14.4935	0.100138
+679.98	14.4896	0.0997902
+680.38	14.4857	0.099473
+680.781	14.4819	0.0991889
+681.181	14.478	0.0989407
+681.582	14.4742	0.0987305
+681.982	14.4705	0.0985566
+682.382	14.4668	0.0984165
+682.783	14.4631	0.0983078
+683.183	14.4594	0.0982277
+683.584	14.4557	0.0981739
+683.984	14.4521	0.0981437
+684.384	14.4485	0.0981345
+684.785	14.4449	0.0981439
+685.185	14.4413	0.0981692
+685.586	14.4378	0.098208
+685.986	14.4342	0.0982576
+686.386	14.4307	0.0983155
+686.787	14.4272	0.0983792
+687.187	14.4236	0.0984461
+687.588	14.4201	0.0985136
+687.988	14.4166	0.0985792
+688.388	14.4131	0.0986404
+688.789	14.4095	0.0986946
+689.189	14.406	0.0987397
+689.59	14.4025	0.0987756
+689.99	14.399	0.098803
+690.39	14.3954	0.0988221
+690.791	14.3919	0.0988336
+691.191	14.3883	0.098838
+691.592	14.3848	0.0988358
+691.992	14.3813	0.0988274
+692.392	14.3777	0.0988133
+692.793	14.3742	0.0987942
+693.193	14.3706	0.0987704
+693.594	14.367	0.0987424
+693.994	14.3635	0.0987109
+694.394	14.3599	0.0986762
+694.795	14.3564	0.0986388
+695.195	14.3528	0.0985993
+695.596	14.3492	0.0985582
+695.996	14.3457	0.098516
+696.396	14.3421	0.0984731

examples/Si.txt

@@ -0,0 +1,300 @@
+#wavelength	Real	Imag	
+101.2	-0.777188	0.470016
+103.3	-0.86132	0.494982
+105.5	-0.953019	0.52318 
+107.8	-1.054431	0.55756
+110.2	-1.160811	0.59274
+112.7	-1.271616	0.63104
+115.3	-1.386816	0.67276
+118.1	-1.505804	0.72072
+121	-1.655375	0.7788
+124	-1.810764	0.84456
+127.2	-2.001376	0.9222
+130.5	-2.169876	1.00264
+134	-2.375296	1.09968
+137.8	-2.620911	1.21844
+141.7	-2.841579	1.34594
+145.9	-3.141004	1.50696
+150.3	-3.412864	1.6872
+155	-3.771516	1.912
+160	-4.142711	2.1714
+165.3	-4.567131	2.48846
+171	-5.000476	2.86752
+177.1	-5.537376	3.3418
+183.7	-6.084864	3.90096
+190.7	-6.735491	4.62462
+198.4	-7.415076	5.59504
+206.6	-8.0704	6.771
+206.6	-7.442181	5.87618
+207.3	-7.499888	6.066816
+208	-7.57302	6.158848
+208.7	-7.489613	6.261684
+209.4	-7.634469	6.34082
+210.1	-7.719435	6.459012
+210.9	-7.738425	6.499712
+211.6	-7.815621	6.62302
+212.3	-7.860344	6.68727
+213	-7.898464	6.76995
+213.8	-8.5748	7.8864
+213.8	-7.988336	6.89997
+214.5	-8.0724	6.972958
+215.3	-8.109304	7.09863
+216	-8.1685	7.184208
+216.8	-8.241779	7.2897
+217.5	-8.292144	7.34432
+218.3	-8.4002	7.49265
+219.1	-8.455979	7.6293
+219.8	-8.549277	7.745036
+220.6	-8.650875	7.8793
+221.4	-9.1289	9.324
+221.4	-8.723427	7.995764
+222.2	-8.819759	8.16684
+223	-8.891625	8.3072
+223.8	-8.985888	8.487666
+224.6	-9.051464	8.66583
+225.4	-9.2799	10.268
+225.4	-9.107604	8.84936
+226.3	-9.160717	9.040956
+227.1	-9.186235	9.261852
+227.9	-9.217444	9.4872
+228.8	-9.194856	9.70751
+229.6	-9.144	11.2832
+229.6	-9.166115	9.902772
+230.5	-9.08742	10.117472
+231.3	-9.014748	10.279136
+232.2	-8.920192	10.414944
+233.1	-8.823808	10.517256
+233.9	-8.6339	11.97
+233.9	-8.749368	10.588774
+234.8	-8.683491	10.60682
+235.7	-8.654236	10.63392
+236.6	-8.651055	10.640608
+237.5	-8.667945	10.656848
+238.4	-8.1212	11.9616
+238.4	-8.724395	10.659012
+239.4	-8.794373	10.667436
+240.3	-8.894652	10.681664
+241.2	-8.99208	10.702462
+242.2	-9.140135	10.721568
+243.1	-8.738	11.7648
+243.1	-9.29	10.773918
+244.1	-9.444501	10.83614
+245	-9.627768	10.910826
+246	-9.819392	10.988544
+247	-10.020328	11.086554
+248	-9.994	12.0288
+248	-10.244325	11.1941
+249	-10.464979	11.3337
+250	-10.695024	11.47712
+251	-10.959844	11.62656
+252	-11.2254	11.802038
+253	-11.504592	11.974306
+254.1	-11.770857	12.185424
+255.1	-12.089301	12.41006
+256.2	-12.400759	12.64104
+257.2	-12.729735	12.907408
+258.3	-13.083477	13.194364
+259.4	-13.450032	13.487726
+260.5	-13.84888	13.833792
+261.6	-14.279832	14.214474
+262.7	-14.715352	14.629014
+263.8	-15.189588	15.092784
+264.9	-15.704064	15.6078
+266.1	-16.236915	16.208012
+267.2	-16.79216	16.888488
+268.4	-17.354771	17.68986
+269.5	-17.93154	18.599728
+270.7	-18.456215	19.618152
+271.9	-18.933201	20.75372
+273.1	-19.343733	22.040644
+274.3	-19.6392	23.441458
+275.5	-19.819323	24.911964
+276.8	-19.88672	26.481312
+278	-19.812436	28.1124
+279.2	-19.61016	29.786162
+280.5	-19.28774	31.548192
+281.8	-18.823936	33.34656
+283.1	-18.216432	35.267074
+284.4	-17.45026	37.291632
+285.7	-16.333269	39.50158
+287	-14.76792	41.869422
+288.3	-12.410629	44.08794
+289.7	-9.455477	45.779436
+291	-6.129339	46.6799
+292.4	-2.931525	46.756908
+293.8	-0.067319	46.24332
+295.2	2.372223	45.350864
+296.6	4.343081	44.27136
+298	5.978304	43.15059
+299.5	7.316385	42.031592
+300.9	8.424748	40.958064
+302.4	9.367959	39.94916
+303.9	10.117216	39.01317
+305.4	10.775596	38.13192
+306.9	11.341924	37.33392
+308.4	11.827725	36.6095
+310	12.233531	35.94174
+311.5	12.63624	35.353522
+313.1	13.010571	34.83954
+314.7	13.33204	34.351722
+316.3	13.655148	33.941264
+317.9	13.96454	33.567072
+319.5	14.254341	33.23902
+321.2	14.565216	32.93995
+322.9	14.891036	32.67936
+324.6	15.212183	32.443944
+326.3	15.529101	32.23366
+328	15.848525	32.038332
+329.7	16.174641	31.87432
+331.5	16.534616	31.74369
+333.3	16.877715	31.615172
+335.1	17.232972	31.534096
+336.9	17.58652	31.477962
+338.8	17.955175	31.442568
+340.6	18.321312	31.459234
+342.5	18.708096	31.51339
+344.4	19.125447	31.638304
+346.3	19.574607	31.834576
+348.3	20.072433	32.125744
+350.2	20.681243	32.532276
+352.2	21.421224	33.07897
+354.2	22.387904	33.81708
+356.3	23.710613	34.696116
+358.3	25.600707	35.635124
+360.4	28.183597	36.314796
+362.5	31.490703	36.346696
+364.7	35.219459	35.28402
+366.8	38.791089	32.88584
+369	41.481055	29.481048
+371.2	43.138341	25.54046
+373.4	43.75836	21.362458
+375.7	43.26564	17.725178
+378	42.130125	14.6187
+380.3	40.735279	12.19428
+382.7	39.227631	10.29508
+385	37.744429	8.83218
+387.5	36.350944	7.63812
+389.9	35.063983	6.673656
+392.4	33.873939	5.90042
+394.9	32.7856	5.238528
+397.4	31.79466	4.704128
+400	30.875131	4.31118
+402.5	30.047024	3.90003
+405.2	29.268159	3.56636
+407.8	28.513832	3.348474
+410.5	27.835975	3.075288
+413.3	27.196923	2.809436
+416.1	26.601871	2.63364
+418.9	26.042345	2.493192
+421.7	25.53138	2.306448
+424.6	25.04556	2.113798
+427.5	24.570312	2.014166
+430.5	24.12942	1.907408
+433.5	23.702159	1.80264
+436.6	23.304336	1.78747
+439.7	22.924781	1.62894
+442.8	22.56444	1.549478
+446	22.237323	1.405964
+449.2	21.898923	1.395236
+452.5	21.586215	1.236368
+455.8	21.281064	1.20913
+459.2	20.986989	1.19158
+462.6	20.712648	1.192886
+466.1	20.430528	1.211896
+469.6	20.190625	1.0788
+473.2	19.930756	1.07184
+476.9	19.723264	0.79956
+480.6	19.49222	0.830208
+484.3	19.273992	0.728906
+488.1	19.064448	0.689986
+492	18.85572	0.668822
+495.9	18.657071	0.63072
+499.9	18.467475	0.627508
+504	18.288373	0.564564
+508.1	18.099841	0.61272
+512.3	17.931625	0.5082
+516.6	17.762625	0.5058
+520.9	17.60328	0.469952
+525.4	17.44452	0.442762
+529.9	17.295432	0.357674
+534.4	17.137575	0.3726
+539.1	16.996825	0.395808
+543.8	16.8573	0.361328
+548.6	16.717985	0.359832
+553.5	16.588305	0.260672
+558.5	16.457805	0.308332
+563.6	16.33674	0.258688
+568.7	16.20752	0.273768
+574	16.095244	0.24072
+579.4	15.97528	0.215838
+584.8	15.863389	0.23898
+590.4	15.752061	0.23814
+596.1	15.649207	0.213624
+601.9	15.546624	0.19715
+607.8	15.452136	0.19655
+613.8	15.350148	0.188064
+619.9	15.256352	0.171864
+626.2	15.154965	0.171292
+632.6	15.069563	0.147516
+639.1	14.976576	0.13932
+645.8	14.883875	0.131172
+652.6	14.799153	0.123104
+659.5	14.722313	0.122784
+666.6	14.638051	0.11478
+673.8	14.554029	0.10682
+681.2	14.477856	0.09893
+688.8	14.409447	0.098696
+696.5	14.3412	0.098462
+704.5	14.27314	0.090672
+712.6	14.197703	0.082896
+720.8	14.145	0.082742
+729.3	14.077404	0.07504
+738	14.024925	0.0749
+746.9	13.957615	0.067248
+756	13.897903	0.067104
+765.3	13.845777	0.059536
+774.9	13.793732	0.059424
+784.7	13.726976	0.05187
+794.8	13.66776	0.051758
+805.1	13.601308	0.044256
+815.7	13.549725	0.044172
+826.6	13.490904	0.03673
+968.6	13.00859475	0.004616631
+984	12.95681319	0.003671549
+999.9	12.9034587	0.002873711
+1016	12.84954523	0.002007391
+1033	12.79274043	0.00128761
+1051	12.73273131	0.000856391
+1069	12.67286325	0.000477026
+1088	12.60982224	0.000298287
+1107	12.54693842	0.000177108
+1120	12.50400321	0.000121642
+1127	12.49039293	9.18886E-05
+1144	12.45737025	5.06231E-05
+1148	12.4522278	4.09337E-05
+1170	12.42396327	1.05743E-05
+1200	12.38547249	1.03532E-05
+1372	12.25490049	9.13077E-06
+1400	12.16335376	8.90722E-06
+1532	12.09926656	7.99332E-06
+1600	12.047841	7.51859E-06
+1696	12.00206736	6.85933E-06
+1800	11.95638084	6.14888E-06
+2000	11.895601	4.79553E-06
+2437.3	11.85700356	1.8676E-06
+2500	11.85011776	1.44849E-06
+2714.4	11.82878449	1.71965E-08
+3000	11.80678321	3.61639E-08
+3303.3	11.78892225	5.62811E-08
+3418.8	11.78342929	6.39375E-08
+3500	11.77931041	6.93173E-08
+3800	11.76819303	8.91925E-08
+4000	11.76078436	3.01787E-07
+5000	11.73816121	1.36359E-06
+5128	11.73649481	1.56219E-06
+5263	11.73473742	2.03481E-06
+5405	11.73288903	2.23332E-06
+5556	11.73092366	2.69893E-06
+5714	11.72886736	2.77404E-06
+5882	11.7266811	2.85597E-06
+6000	11.72514564	2.91307E-06

examples/coating-flow.py

@@ -94,28 +94,55 @@ if __name__ == '__main__':
             for i in range(len(nvalues)):
                 r += [r[i] + tl]
 
-            x = np.ones((1, len(nvalues) + 1), dtype = np.float64)
-            m = np.ones((1, len(nvalues) + 1), dtype = np.complex128)
+            x = np.ones((len(nvalues) + 1), dtype = np.float64)
+            m = np.ones((len(nvalues) + 1), dtype = np.complex128)
 
-            x[0] = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
-            m[0] = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
+            x = 2.0*np.pi*np.array(r, dtype = np.float64)/wl
+            m = np.array([ms] + nvalues[:, 1].tolist(), dtype = np.complex128)
             print(x,m)
 
-            factor = 2
+            factor = 2.91*x[0]/x[-1]
+            print factor
             comment='PEC-'+basename
-            WL_units='cm'
+            WL_units=''
             #flow_total = 39
-            flow_total = 25
-            crossplane='XZ'
+            # flow_total = 23 #SV False
+            flow_total = 24
+            #flow_total = 4
+            #crossplane='XZ'
+            crossplane='XYZ'
             #crossplane='YZ'
             #crossplane='XY'
 
             # Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
-            field_to_plot='Pabs'
+            #field_to_plot='Pabs'
             #field_to_plot='Eabs'
-            #field_to_plot='angleEx'
-            fieldplot(x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts,
-                      factor, flow_total, pl=0, outline_width=2.0)
+            
+            field_to_plot='angleEx'
+            #field_to_plot='angleHy'
+            print "x =", x
+            print "m =", m
+
+            import matplotlib.pyplot as plt
+            plt.rcParams.update({'font.size': 16})
+            fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+            fig.tight_layout()
+            fieldplot(fig, axs, x,m, wl, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+                      subplot_label=' ',is_flow_extend=False
+                      , outline_width=1.5
+                      , pl=0 #PEC layer starts the design
+                      )
+            # fieldplot(fig, axs, x[0],m[0], wl, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+            #           subplot_label=' ' ,is_flow_extend=False
+            #           , outline_width=1.5
+            #           , pl=0 #PEC layer starts the design
+            #           )
+            fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+            plt.savefig(comment+"-R"+str(int(round(x[-1]*wl/2.0/np.pi)))+"-"+crossplane+"-"
+                        +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+            plt.draw()
+            plt.clf()
+            plt.close()
 
 
         print "Done!!"

examples/coating-flow.sh

@@ -5,4 +5,4 @@
 #python coating-flow.py -w 3.75 -t 0.62 -f index-sv.dat -n 1501
 # python coating-flow.py -w 3.75 -t 2.40 -f index-ch.dat -n 501
 
-python coating-flow.py -w 0.532 -t 0.128 -f index-in-glass.dat -n 1501
+python coating-flow.py -w 0.532 -t 0.128 -f index-in-glass.dat -n 151

examples/example-get-Mie.cc

@@ -0,0 +1,202 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program is free software: you can redistribute it and/or modify          //
+//    it under the terms of the GNU General Public License as published by          //
+//    the Free Software Foundation, either version 3 of the License, or             //
+//    (at your option) any later version.                                           //
+//                                                                                  //
+//    This program is distributed in the hope that it will be useful,               //
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
+//    GNU General Public License for more details.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//                                                                                  //
+//    You should have received a copy of the GNU General Public License             //
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+//   This program returns expansion coefficents of Mie series
+#include <complex>
+#include <cstdio>
+#include <string>
+#include "../src/nmie-applied.h"
+#include "./read-spectra.h"
+template<class T> inline T pow2(const T value) {return value*value;}
+int main(int argc, char *argv[]) {
+  try {
+    read_spectra::ReadSpectra Si_index, Ag_index;
+    read_spectra::ReadSpectra plot_core_index_, plot_TiN_;
+    std::string core_filename("Si-int.txt");
+    //std::string core_filename("Ag.txt");
+    //std::string TiN_filename("TiN.txt");
+    std::string TiN_filename("Ag-int.txt");
+    //std::string TiN_filename("Si.txt");
+    std::string shell_filename(core_filename);
+
+    nmie::MultiLayerMieApplied multi_layer_mie;  
+    const std::complex<double> epsilon_Si(18.4631066585, 0.6259727805);
+    const std::complex<double> epsilon_Ag(-8.5014154589, 0.7585845411);
+    const std::complex<double> index_Si = std::sqrt(epsilon_Si);
+    const std::complex<double> index_Ag = std::sqrt(epsilon_Ag);
+    double WL=500; //nm
+    double core_width = 5.27; //nm Si
+    double inner_width = 8.22; //nm Ag
+    double outer_width = 67.91; //nm  Si
+    bool isSiAgSi = true;
+    double delta_width = 25.0;
+    //bool isSiAgSi = false;
+    if (isSiAgSi) {
+      core_width = 5.27; //nm Si
+      inner_width = 8.22; //nm Ag
+      outer_width = 67.91; //nm  Si
+      multi_layer_mie.AddTargetLayer(core_width, index_Si);
+      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+    } else {
+      inner_width = 31.93; //nm Ag
+      outer_width = 4.06; //nm  Si
+      multi_layer_mie.AddTargetLayer(inner_width, index_Ag);
+      multi_layer_mie.AddTargetLayer(outer_width+delta_width, index_Si);
+    }
+
+    for (int dd = 0; dd<50; ++dd) {
+      delta_width = dd;
+    FILE *fp;
+    std::string fname = "absorb-layered-spectra-d"+std::to_string(dd)+".dat";
+    fp = fopen(fname.c_str(), "w");
+
+    multi_layer_mie.SetWavelength(WL);
+    multi_layer_mie.RunMieCalculation();
+    double Qabs = multi_layer_mie.GetQabs();
+    printf("Qabs = %g\n", Qabs);
+    std::vector< std::vector<std::complex<double> > > aln, bln, cln, dln;
+    multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+    std::string str = std::string("#WL ");
+    for (int l = 0; l<aln.size(); ++l) {
+      for (int n = 0; n<3; ++n) {
+	str+="|a|^2+|d|^2_ln"+std::to_string(l)+std::to_string(n)+" "
+	  + "|b|^2+|c|^2_ln"+std::to_string(l)+std::to_string(n)+" ";
+      }
+    }
+    str+="\n";
+    fprintf(fp, "%s", str.c_str());
+    fprintf(fp, "# |a|+|d|");
+    str.clear();
+
+    double from_WL = 400;
+    double to_WL = 600;
+    int total_points = 401;
+    double delta_WL = std::abs(to_WL - from_WL)/(total_points-1.0);
+    Si_index.ReadFromFile(core_filename).ResizeToComplex(from_WL, to_WL, total_points)
+      .ToIndex();
+    Ag_index.ReadFromFile(TiN_filename).ResizeToComplex(from_WL, to_WL, total_points)
+      .ToIndex();
+    auto Si_data = Si_index.GetIndex();
+    auto Ag_data = Ag_index.GetIndex();
+    for (int i=0; i < Si_data.size(); ++i) {
+      const double& WL = Si_data[i].first;
+      const std::complex<double>& Si = Si_data[i].second;
+      const std::complex<double>& Ag = Ag_data[i].second;
+      str+=std::to_string(WL);
+      multi_layer_mie.ClearTarget();      
+      if (isSiAgSi) {
+	multi_layer_mie.AddTargetLayer(core_width, Si);
+	multi_layer_mie.AddTargetLayer(inner_width, Ag);
+	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
+      } else {
+	inner_width = 31.93; //nm Ag
+	outer_width = 4.06; //nm  Si
+	multi_layer_mie.AddTargetLayer(inner_width, Ag);
+	multi_layer_mie.AddTargetLayer(outer_width+delta_width, Si);
+      }
+      multi_layer_mie.SetWavelength(WL);
+      multi_layer_mie.RunMieCalculation();
+      multi_layer_mie.GetQabs();
+      multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+      for (int l = 0; l<aln.size(); ++l) {
+	for (int n = 0; n<3; ++n) {
+	  str+=" "+std::to_string(pow2(std::abs(aln[l][n]))+
+				   pow2(std::abs(dln[l][n])))
+	    + " "
+	    + std::to_string(pow2(std::abs(bln[l][n]))
+				 + pow2(std::abs(cln[l][n])) );
+
+	  // str+=" "+std::to_string(aln[l][n].real() - pow2(std::abs(aln[l][n]))
+	  // 			  +dln[l][n].real() - pow2(std::abs(dln[l][n])))
+	  //   + " "
+	  //   + std::to_string(bln[l][n].real() - pow2(std::abs(bln[l][n]))
+	  // 			  +cln[l][n].real() - pow2(std::abs(cln[l][n])) );
+	}
+      }
+      str+="\n";
+      fprintf(fp, "%s", str.c_str());
+      str.clear();
+    }
+
+    fclose(fp);
+    }
+
+    // WL = 500;
+    // multi_layer_mie.SetWavelength(WL);
+    // multi_layer_mie.RunMieCalculation();
+    // multi_layer_mie.GetQabs();
+    // multi_layer_mie.GetExpanCoeffs(aln, bln, cln, dln);
+
+    // printf("\n Scattering");
+    // for (int l = 0; l<aln.size(); ++l) {
+    //   int n = 0;
+    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    //   n = 1;
+    //   printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+    //   printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+    //   printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+    //   printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    //   // n = 2;
+    //   // printf("aln[%i][%i] = %g, %gi\n", l, n+1, aln[l][n].real(), aln[l][n].imag());
+    //   // printf("bln[%i][%i] = %g, %gi\n", l, n+1, bln[l][n].real(), bln[l][n].imag());
+    //   // printf("cln[%i][%i] = %g, %gi\n", l, n+1, cln[l][n].real(), cln[l][n].imag());
+    //   // printf("dln[%i][%i] = %g, %gi\n", l, n+1, dln[l][n].real(), dln[l][n].imag());
+    // }
+    // printf("\n Absorbtion\n");
+    // for (int l = 0; l<aln.size(); ++l) {
+    //   if (l == aln.size()-1) printf(" Total ");
+    //   printf("===== l=%i   =====\n", l);
+    //   int n = 0;
+    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
+    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
+    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
+    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
+    //   n = 1;
+    //   printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
+    //   printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
+    //   printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
+    //   printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
+    //   // n = 2;
+    //   // printf("aln[%i][%i] = %g\n", l, n+1, aln[l][n].real() - pow2(std::abs(aln[l][n])));
+    //   // printf("bln[%i][%i] = %g\n", l, n+1, bln[l][n].real() - pow2(std::abs(bln[l][n])));
+    //   // printf("cln[%i][%i] = %g\n", l, n+1, cln[l][n].real() - pow2(std::abs(cln[l][n])));
+    //   // printf("dln[%i][%i] = %g\n", l, n+1, dln[l][n].real() - pow2(std::abs(dln[l][n])));
+    // }
+
+
+  } catch( const std::invalid_argument& ia ) {
+    // Will catch if  multi_layer_mie fails or other errors.
+    std::cerr << "Invalid argument: " << ia.what() << std::endl;
+    return -1;
+  }  
+    return 0;
+}
+
+

examples/field-Ag-flow.py

@@ -87,5 +87,27 @@ crossplane='XZ'
 field_to_plot='Pabs'
 #field_to_plot='angleEx'
 
-fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
 
+import matplotlib.pyplot as plt
+fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+fig.tight_layout()
+fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+          subplot_label=' ',is_flow_extend=False)
+
+#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+
+# for ax in axs:
+#     ax.locator_params(axis='x',nbins=5)
+#     ax.locator_params(axis='y',nbins=5)
+
+fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+
+plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
+                    +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+
+plt.draw()
+
+#    plt.show()
+
+plt.clf()
+plt.close()

examples/field-SiAgSi-flow.py

@@ -2,6 +2,7 @@
 # -*- coding: UTF-8 -*-
 #
 #    Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
+#    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>
 #
 #    This file is part of python-scattnlay
 #
@@ -33,43 +34,11 @@
 import scattnlay
 from scattnlay import fieldnlay
 from scattnlay import scattnlay
+from fieldplot import fieldplot
 import numpy as np
 import cmath
 
 
-def get_index(array,value):
-    idx = (np.abs(array-value)).argmin()
-    return idx
-
-#Ec = np.resize(Ec, (npts, npts)).T
-
-
-def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
-    # Initial position
-    flow_x = [a]
-    flow_z = [b]
-    for n in range(0, nmax):
-        #Get the next position
-        #1. Find Poynting vector and normalize it
-        x_pos = flow_x[-1]
-        z_pos = flow_z[-1]
-        x_idx = get_index(scale_x, x_pos)
-        z_idx = get_index(scale_z, z_pos)
-        S=np.cross(Ec[npts*z_idx+x_idx], np.conjugate(Hc[npts*z_idx+x_idx]) ).real
-        Snorm=S/np.linalg.norm(S)
-        #2. Evaluate displacement = half of the discrete and new position
-        dpos = abs(scale_z[0]-scale_z[1])/4.0
-        dx = dpos*Snorm[0]
-        dz = dpos*Snorm[2]
-        x_pos = x_pos+dx
-        z_pos = z_pos+dz
-        #3. Save result
-        flow_x.append(x_pos)
-        flow_z.append(z_pos)
-    return flow_x, flow_z
-
-
-
 epsilon_Si = 13.64 + 0.047j
 epsilon_Ag = -28.05 + 1.525j
 
@@ -103,135 +72,56 @@ outer_r = inner_r+outer_width
 # n2 = 0.565838 + 7.23262j
 nm = 1.0
 
-x = np.ones((1, 3), dtype = np.float64)
-x[0, 0] = 2.0*np.pi*core_r/WL
-x[0, 1] = 2.0*np.pi*inner_r/WL
-x[0, 2] = 2.0*np.pi*outer_r/WL
+x = np.ones((3), dtype = np.float64)
+x[0] = 2.0*np.pi*core_r/WL
+x[1] = 2.0*np.pi*inner_r/WL
+x[2] = 2.0*np.pi*outer_r/WL
 
-m = np.ones((1, 3), dtype = np.complex128)
-m[0, 0] = index_Si/nm
-m[0, 1] = index_Ag/nm
-m[0, 2] = index_Si/nm
+m = np.ones((3), dtype = np.complex128)
+m[0] = index_Si/nm
+m[1] = index_Ag/nm
+m[2] = index_Si/nm
 
 print "x =", x
 print "m =", m
 
-npts = 241
-
+npts = 501
 factor=2.2
-scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
-
-coordX, coordZ = np.meshgrid(scan, scan)
-coordX.resize(npts*npts)
-coordZ.resize(npts*npts)
-coordY = np.zeros(npts*npts, dtype = np.float64)
-
-coord = np.vstack((coordX, coordY, coordZ)).transpose()
-
-terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
-terms, E, H = fieldnlay(x, m, coord)
-Er = np.absolute(E)
-Hr = np.absolute(H)
-
-# |E|/|Eo|
-Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
-Ec = E[0, :, :]
-Hc = H[0, :, :]
-Eangle = np.angle(E[0, :, 0])/np.pi*180
-
-P=[]
-for n in range(0, len(E[0])):
-    P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
-
-Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
-Hangle = np.angle(H[0, :, 1])/np.pi*180
-
-
-
-try:
-    import matplotlib.pyplot as plt
-    from matplotlib import cm
-    from matplotlib.colors import LogNorm
-
-    min_tick = 0.0
-    max_tick = 1.0
-
-    # Eabs_data = np.resize(P, (npts, npts)).T
-    Eabs_data = np.resize(Eabs, (npts, npts)).T
-    # Eangle_data = np.resize(Eangle, (npts, npts)).T
-    # Habs_data = np.resize(Habs, (npts, npts)).T
-    # Hangle_data = np.resize(Hangle, (npts, npts)).T
-
-    fig, ax = plt.subplots(1,1)#, sharey=True, sharex=True)
-    #fig.tight_layout()
-    # Rescale to better show the axes
-    scale_x = np.linspace(min(coordX)*WL/2.0/np.pi/nm, max(coordX)*WL/2.0/np.pi/nm, npts)
-    scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi/nm, max(coordZ)*WL/2.0/np.pi/nm, npts)
-
-    # Define scale ticks
-    min_tick = min(min_tick, np.amin(Eabs_data))
-    max_tick = max(max_tick, np.amax(Eabs_data))
-    #max_tick = 5
-    # scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
-    scale_ticks = np.linspace(min_tick, max_tick, 11)
-
-    # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
-    ax.set_title('Eabs')
-    cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = cm.jet,
-                    origin = 'lower'
-                    , vmin = min_tick, vmax = max_tick
-                    , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
-                    #,norm = LogNorm()
-                    )
-    ax.axis("image")
-
-    # Add colorbar
-    cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
-    cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
-    pos = list(cbar.ax.get_position().bounds)
-    fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
-
-    plt.xlabel('Z, nm')
-    plt.ylabel('X, nm')
-
-    # This part draws the nanoshell
-    from matplotlib import patches
-    s1 = patches.Arc((0, 0), 2.0*core_r, 2.0*core_r,  angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    s2 = patches.Arc((0, 0), 2.0*inner_r, 2.0*inner_r, angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    s3 = patches.Arc((0, 0), 2.0*outer_r, 2.0*outer_r, angle=0.0, zorder=2,
-                     theta1=0.0, theta2=360.0, linewidth=1, color='black')
-    ax.add_patch(s1)
-    ax.add_patch(s2) 
-    ax.add_patch(s3) 
-
-    from matplotlib.path import Path
-    #import matplotlib.patches as patches
-
-    flow_total = 21
-    for flow in range(0,flow_total):
-        flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
-                                 min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
-                                                    min(scale_z), npts*12)
-        verts = np.vstack((flow_z, flow_x)).transpose().tolist()
-        codes = [Path.CURVE4]*len(verts)
-        #codes = [Path.LINETO]*len(verts)
-        codes[0] = Path.MOVETO
-        path = Path(verts, codes)
-        patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
-        ax.add_patch(patch)
-
- 
-    plt.savefig("SiAgSi-flow.png")
-    plt.draw()
-
-    plt.show()
-
-    plt.clf()
-    plt.close()
-finally:
-    print("Qabs = "+str(Qabs));
-#
+flow_total = 21
+
+
+crossplane='XZ'
+#crossplane='YZ'
+#crossplane='XY'
+
+# Options to plot: Eabs, Habs, Pabs, angleEx, angleHy
+field_to_plot='Eabs'
+#field_to_plot='angleEx'
+comment='SiAgSi-absorber-flow'
+WL_units='nm'
+
+import matplotlib.pyplot as plt
+fig, axs = plt.subplots(1,1)#, sharey=True, sharex=True)
+fig.tight_layout()
+fieldplot(fig, axs, x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total,
+          subplot_label=' ',is_flow_extend=False, outline_width=1.5)
+
+#fieldplot(x,m, WL, comment, WL_units, crossplane, field_to_plot, npts, factor, flow_total, is_flow_extend=False)
+
+# for ax in axs:
+#     ax.locator_params(axis='x',nbins=5)
+#     ax.locator_params(axis='y',nbins=5)
+
+fig.subplots_adjust(hspace=0.3, wspace=-0.1)
+
+plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
+                    +field_to_plot+".pdf",pad_inches=0.02, bbox_inches='tight')
+
+plt.draw()
+
+#    plt.show()
+
+plt.clf()
+plt.close()
 
 

examples/fieldplot.py

@@ -39,6 +39,7 @@ def unit_vector(vector):
     """ Returns the unit vector of the vector.  """
     return vector / np.linalg.norm(vector)
 
+
 def angle_between(v1, v2):
     """ Returns the angle in radians between vectors 'v1' and 'v2'::
 
@@ -59,66 +60,70 @@ def angle_between(v1, v2):
             return np.pi
     return angle
 ###############################################################################
+
+
 def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
     # Initial position
     flow_x = [x0]
     flow_y = [y0]
     flow_z = [z0]
-    max_step = x[-1]/3
-    min_step = x[0]/2000
+    max_step = x[-1] / 3
+    min_step = x[0] / 2000
 #    max_step = min_step
-    step = min_step*2.0
+    step = min_step * 2.0
     if max_step < min_step:
         max_step = min_step
     coord = np.vstack(([flow_x[-1]], [flow_y[-1]], [flow_z[-1]])).transpose()
-    terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord,pl=pl)
+    terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord, pl=pl)
     Ec, Hc = E[0, 0, :], H[0, 0, :]
     S = np.cross(Ec, Hc.conjugate()).real
-    Snorm_prev = S/np.linalg.norm(S)
+    Snorm_prev = S / np.linalg.norm(S)
     Sprev = S
     length = 0
     dpos = step
     count = 0
     while length < max_length:
         count = count + 1
-        if (count>3000): # Limit length of the absorbed power streamlines
+        if (count > 4000):  # Limit length of the absorbed power streamlines
             break
-        if step<max_step:
-                step = step*2.0
+        if step < max_step:
+            step = step * 2.0
         r = np.sqrt(flow_x[-1]**2 + flow_y[-1]**2 + flow_z[-1]**2)
         while step > min_step:
-            #Evaluate displacement from previous poynting vector
+            # Evaluate displacement from previous poynting vector
             dpos = step
-            dx = dpos*Snorm_prev[0];
-            dy = dpos*Snorm_prev[1];
-            dz = dpos*Snorm_prev[2];
-            #Test the next position not to turn\chang size for more than max_angle
-            coord = np.vstack(([flow_x[-1]+dx], [flow_y[-1]+dy], [flow_z[-1]+dz])).transpose()
-            terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord,pl=pl)
+            dx = dpos * Snorm_prev[0]
+            dy = dpos * Snorm_prev[1]
+            dz = dpos * Snorm_prev[2]
+            # Test the next position not to turn\chang size for more than
+            # max_angle
+            coord = np.vstack(([flow_x[-1] + dx], [flow_y[-1] + dy],
+                               [flow_z[-1] + dz])).transpose()
+            terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord, pl=pl)
             Ec, Hc = E[0, 0, :], H[0, 0, :]
-            Eth = max(np.absolute(Ec))/1e10
-            Hth = max(np.absolute(Hc))/1e10
-            for i in xrange(0,len(Ec)):
+            Eth = max(np.absolute(Ec)) / 1e10
+            Hth = max(np.absolute(Hc)) / 1e10
+            for i in xrange(0, len(Ec)):
                 if abs(Ec[i]) < Eth:
-                    Ec[i] = 0+0j
+                    Ec[i] = 0 + 0j
                 if abs(Hc[i]) < Hth:
-                    Hc[i] = 0+0j
+                    Hc[i] = 0 + 0j
             S = np.cross(Ec, Hc.conjugate()).real
             if not np.isfinite(S).all():
                 break
-            Snorm = S/np.linalg.norm(S)
-            diff = (S-Sprev)/max(np.linalg.norm(S), np.linalg.norm(Sprev))
-            if np.linalg.norm(diff)<max_angle:
-            # angle = angle_between(Snorm, Snorm_prev)
-            # if abs(angle) < max_angle:
+            Snorm = S / np.linalg.norm(S)
+            diff = (S - Sprev) / max(np.linalg.norm(S), np.linalg.norm(Sprev))
+            if np.linalg.norm(diff) < max_angle:
+                # angle = angle_between(Snorm, Snorm_prev)
+                # if abs(angle) < max_angle:
                 break
-            step = step/2.0
-        #3. Save result
+            step = step / 2.0
+        # 3. Save result
         Sprev = S
         Snorm_prev = Snorm
-        dx = dpos*Snorm_prev[0];
-        dy = dpos*Snorm_prev[1];
-        dz = dpos*Snorm_prev[2];
+        dx = dpos * Snorm_prev[0]
+        dy = dpos * Snorm_prev[1]
+        dz = dpos * Snorm_prev[2]
         length = length + step
         flow_x.append(flow_x[-1] + dx)
         flow_y.append(flow_y[-1] + dy)
@@ -129,7 +134,7 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m, pl):
 ###############################################################################
 def GetField(crossplane, npts, factor, x, m, pl):
     """
-    crossplane: XZ, YZ, XY
+    crossplane: XZ, YZ, XY, or XYZ (half is XZ, half is YZ)
     npts: number of point in each direction
     factor: ratio of plotting size to outer size of the particle
     x: size parameters for particle layers
@@ -140,172 +145,226 @@ def GetField(crossplane, npts, factor, x, m, pl):
 
     if crossplane=='XZ':
         coordX, coordZ = np.meshgrid(scan, scan)
-        coordX.resize(npts*npts)
-        coordZ.resize(npts*npts)
+        coordX.resize(npts * npts)
+        coordZ.resize(npts * npts)
         coordY = zero
         coordPlot1 = coordX
         coordPlot2 = coordZ
-    elif crossplane=='YZ':
+    elif crossplane == 'YZ':
         coordY, coordZ = np.meshgrid(scan, scan)
-        coordY.resize(npts*npts)
-        coordZ.resize(npts*npts)
+        coordY.resize(npts * npts)
+        coordZ.resize(npts * npts)
         coordX = zero
         coordPlot1 = coordY
         coordPlot2 = coordZ
-    elif crossplane=='XY':
+    elif crossplane == 'XY':
         coordX, coordY = np.meshgrid(scan, scan)
-        coordX.resize(npts*npts)
-        coordY.resize(npts*npts)
+        coordX.resize(npts * npts)
+        coordY.resize(npts * npts)
         coordZ = zero
         coordPlot1 = coordY
         coordPlot2 = coordX
-        
+    elif crossplane=='XYZ':
+        coordX, coordZ = np.meshgrid(scan, scan)
+        coordY, coordZ = np.meshgrid(scan, scan)
+        coordPlot1, coordPlot2 = np.meshgrid(scan, scan)
+        coordPlot1.resize(npts * npts)
+        coordPlot2.resize(npts * npts)
+        half=npts//2
+        # coordX = np.copy(coordX)
+        # coordY = np.copy(coordY)
+        coordX[:,:half]=0
+        coordY[:,half:]=0
+        coordX.resize(npts*npts)
+        coordY.resize(npts*npts)
+        coordZ.resize(npts*npts)
+    else:
+        print("Unknown crossplane")
+        import sys
+        sys.exit()
+
     coord = np.vstack((coordX, coordY, coordZ)).transpose()
     terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord, pl=pl)
     Ec = E[0, :, :]
     Hc = H[0, :, :]
-    P=[]
-    P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], np.conjugate(Hc[n]))).real, range(0, len(E[0]))))
+    P = []
+    P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real,
+                     range(0, len(E[0]))))
 
     # for n in range(0, len(E[0])):
     #     P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
     return Ec, Hc, P, coordPlot1, coordPlot2
 ###############################################################################
-def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot='Pabs',npts=101, factor=2.1, flow_total=11, is_flow_extend=True, pl=-1, outline_width=1):
+
+
+def fieldplot(fig, ax, x, m, WL, comment='', WL_units=' ', crossplane='XZ',
+              field_to_plot='Pabs', npts=101, factor=2.1, flow_total=11,
+              is_flow_extend=True, pl=-1, outline_width=1, subplot_label=' '):
     Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m, pl)
     Er = np.absolute(Ec)
     Hr = np.absolute(Hc)
     try:
-        import matplotlib.pyplot as plt
         from matplotlib import cm
         from matplotlib.colors import LogNorm
+
         if field_to_plot == 'Pabs':
-            Eabs_data = np.resize(P, (npts, npts)).T 
-            label = r'$\operatorname{Re}(E \times H^*)$'
+            Eabs_data = np.resize(P, (npts, npts)).T
+            label = r'$\operatorname{Re}(E \times H)$'
         elif field_to_plot == 'Eabs':
-            Eabs = np.sqrt(Er[ :, 0]**2 + Er[ :, 1]**2 + Er[ :, 2]**2)
-            Eabs_data = np.resize(Eabs, (npts, npts)).T 
+            Eabs = np.sqrt(Er[:, 0]**2 + Er[:, 1]**2 + Er[:, 2]**2)
+            Eabs_data = np.resize(Eabs, (npts, npts)).T
             label = r'$|E|$'
         elif field_to_plot == 'Habs':
-            Habs= np.sqrt(Hr[ :, 0]**2 + Hr[ :, 1]**2 + Hr[ :, 2]**2)
-            Eabs_data = np.resize(Habs, (npts, npts)).T 
+            Habs = np.sqrt(Hr[:, 0]**2 + Hr[:, 1]**2 + Hr[:, 2]**2)
+            Eabs_data = np.resize(Habs, (npts, npts)).T
             label = r'$|H|$'
         elif field_to_plot == 'angleEx':
-            Eangle = np.angle(Ec[ :, 0])/np.pi*180
-            Eabs_data = np.resize(Eangle, (npts, npts)).T 
+            Eangle = np.angle(Ec[:, 0]) / np.pi * 180
+            Eabs_data = np.resize(Eangle, (npts, npts)).T
             label = r'$arg(E_x)$'
         elif field_to_plot == 'angleHy':
-            Hangle = np.angle(Hc[ :, 1])/np.pi*180
-            Eabs_data = np.resize(Hangle, (npts, npts)).T 
+            Hangle = np.angle(Hc[:, 1]) / np.pi * 180
+            Eabs_data = np.resize(Hangle, (npts, npts)).T
             label = r'$arg(H_y)$'
 
-        fig, ax = plt.subplots(1,1)
         # Rescale to better show the axes
-        scale_x = np.linspace(min(coordX)*WL/2.0/np.pi, max(coordX)*WL/2.0/np.pi, npts)
-        scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi, max(coordZ)*WL/2.0/np.pi, npts)
+        scale_x = np.linspace(
+            min(coordX) * WL / 2.0 / np.pi, max(coordX) * WL / 2.0 / np.pi, npts)
+        scale_z = np.linspace(
+            min(coordZ) * WL / 2.0 / np.pi, max(coordZ) * WL / 2.0 / np.pi, npts)
 
         # Define scale ticks
         min_tick = np.amin(Eabs_data[~np.isnan(Eabs_data)])
+        #min_tick = 0.1
         max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
-        scale_ticks = np.linspace(min_tick, max_tick, 6)
-
+        #max_tick = 60
+        scale_ticks = np.linspace(min_tick, max_tick, 5)
+        #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
+        #scale_ticks = [0.1,0.3,1,3,10, max_tick]
         # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
         ax.set_title(label)
-        my_cmap = cm.jet
-        if not (field_to_plot == 'angleEx' or field_to_plot == 'angleHy'):
-            my_cmap.set_under('w')
-        cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = my_cmap,
-                        origin = 'lower'
-                        , vmin = min_tick+max_tick*1e-15, vmax = max_tick
-                        , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
-                        #,norm = LogNorm()
+        # build a rectangle in axes coords
+        ax.annotate(subplot_label, xy=(0.0, 1.1), xycoords='axes fraction',  # fontsize=10,
+                    horizontalalignment='left', verticalalignment='top')
+        # ax.text(right, top, subplot_label,
+        #         horizontalalignment='right',
+        #         verticalalignment='bottom',
+        #         transform=ax.transAxes)
+        cax = ax.imshow(Eabs_data, interpolation='nearest', cmap=cm.jet,
+                        origin='lower', vmin=min_tick, vmax=max_tick, extent=(min(scale_x), max(scale_x), min(scale_z), max(scale_z))
+                        # ,norm = LogNorm()
                         )
         ax.axis("image")
 
         # Add colorbar
-        cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
-        cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
-        pos = list(cbar.ax.get_position().bounds)
+        cbar = fig.colorbar(cax, ticks=[a for a in scale_ticks], ax=ax)
+        # vertically oriented colorbar
+        if 'angle' in field_to_plot:
+            cbar.ax.set_yticklabels(['%3.0f' % (a) for a in scale_ticks])
+        else:
+            cbar.ax.set_yticklabels(['%3.2f' % (a) for a in scale_ticks])
+        # pos = list(cbar.ax.get_position().bounds)
         #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
-        if crossplane=='XZ':
-            plt.xlabel('Z, '+WL_units)
-            plt.ylabel('X, '+WL_units)
-        elif crossplane=='YZ':
-            plt.xlabel('Z, '+WL_units)
-            plt.ylabel('Y, '+WL_units)
-        elif crossplane=='XY':
-            plt.xlabel('Y, '+WL_units)
-            plt.ylabel('X, '+WL_units)
-
-
+        lp2 = -10.0
+        lp1 = -1.0
+        if crossplane == 'XZ':
+            ax.set_xlabel('Z, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('X, ' + WL_units, labelpad=lp2)
+        elif crossplane == 'YZ':
+            ax.set_xlabel('Z, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('Y, ' + WL_units, labelpad=lp2)
+        elif crossplane=='XYZ':
+            ax.set_xlabel(r'$Z,\lambda$'+WL_units)
+            ax.set_ylabel(r'$Y:X,\lambda$'+WL_units)
+        elif crossplane == 'XY':
+            ax.set_xlabel('Y, ' + WL_units, labelpad=lp1)
+            ax.set_ylabel('X, ' + WL_units, labelpad=lp2)
         # # This part draws the nanoshell
         from matplotlib import patches
         from matplotlib.path import Path
-        x_edge = (x[-1], x[0])
-        for xx in x_edge:
-            r= xx*WL/2.0/np.pi
-            s1 = patches.Arc((0, 0), 2.0*r, 2.0*r,  angle=0.0, zorder=1.8,
+        for xx in x:
+            r = xx * WL / 2.0 / np.pi
+            s1 = patches.Arc((0, 0), 2.0 * r, 2.0 * r,  angle=0.0, zorder=1.8,
                              theta1=0.0, theta2=360.0, linewidth=outline_width, color='black')
             ax.add_patch(s1)
-        if (crossplane=='XZ' or crossplane=='YZ') and flow_total>0:
+        #
+        # for flow in range(0,flow_total):
+        #     flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
+        #                              min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
+        #                              min(scale_z),
+        #                              #0.0,
+        #                              npts*16)
+        #     verts = np.vstack((flow_z, flow_x)).transpose().tolist()
+        #     #codes = [Path.CURVE4]*len(verts)
+        #     codes = [Path.LINETO]*len(verts)
+        #     codes[0] = Path.MOVETO
+        #     path = Path(verts, codes)
+        #     patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='yellow')
+        #     ax.add_patch(patch)
+        if (not crossplane == 'XY') and flow_total > 0:
 
             from matplotlib.path import Path
-            scanSP = np.linspace(-factor*x[-1], factor*x[-1], npts)
-            min_SP = -factor*x[-1]
-            step_SP = 2.0*factor*x[-1]/(flow_total-1)
+            scanSP = np.linspace(-factor * x[-1], factor * x[-1], npts)
+            min_SP = -factor * x[-1]
+            step_SP = 2.0 * factor * x[-1] / (flow_total - 1)
             x0, y0, z0 = 0, 0, 0
-            max_length=factor*x[-1]*8
-            #max_length=factor*x[-1]*4
-            max_angle = np.pi/160
+            max_length = factor * x[-1] * 10
+            # max_length=factor*x[-1]*5
+            max_angle = np.pi / 160
             if is_flow_extend:
-                rg = range(0,flow_total*2+1)
+                rg = range(0, flow_total * 5 + 1)
             else:
-                rg = range(0,flow_total)
+                rg = range(0, flow_total)
             for flow in rg:
+                if is_flow_extend:
+                    f = min_SP*2 + flow*step_SP
+                else:
+                    f = min_SP + flow*step_SP
                 if crossplane=='XZ':
-                    if is_flow_extend:
-                        x0 = min_SP*2 + flow*step_SP
-                    else:
-                        x0 = min_SP + flow*step_SP
-                    z0 = min_SP
-                    #y0 = x[-1]/20 
+                    x0 = f
                 elif crossplane=='YZ':
-                    if is_flow_extend:
-                        y0 = min_SP*2 + flow*step_SP
+                    y0 = f
+                elif crossplane=='XYZ':
+                    x0 = 0
+                    y0 = 0
+                    if f > 0:
+                        x0 = f
                     else:
-                        y0 = min_SP + flow*step_SP
-                    z0 = min_SP
-                    #x0 = x[-1]/20
-                flow_xSP, flow_ySP, flow_zSP = GetFlow3D(x0, y0, z0, max_length, max_angle, x, m,pl)
-                if crossplane=='XZ':
-                    flow_z_plot = flow_zSP*WL/2.0/np.pi
-                    flow_f_plot = flow_xSP*WL/2.0/np.pi
-                elif crossplane=='YZ':
-                    flow_z_plot = flow_zSP*WL/2.0/np.pi
-                    flow_f_plot = flow_ySP*WL/2.0/np.pi
+                        y0 = f
+                z0 = min_SP
+                    # x0 = x[-1]/20
+                flow_xSP, flow_ySP, flow_zSP = GetFlow3D(
+                    x0, y0, z0, max_length, max_angle, x, m, pl)
+                if crossplane == 'XZ':
+                    flow_z_plot = flow_zSP * WL / 2.0 / np.pi
+                    flow_f_plot = flow_xSP * WL / 2.0 / np.pi
+                elif crossplane == 'YZ':
+                    flow_z_plot = flow_zSP * WL / 2.0 / np.pi
+                    flow_f_plot = flow_ySP * WL / 2.0 / np.pi
+                elif crossplane=='XYZ':
+                    if f > 0:
+                        flow_z_plot = flow_zSP*WL/2.0/np.pi
+                        flow_f_plot = flow_xSP*WL/2.0/np.pi
+                    else:
+                        flow_z_plot = flow_zSP*WL/2.0/np.pi
+                        flow_f_plot = flow_ySP*WL/2.0/np.pi
 
-                verts = np.vstack((flow_z_plot, flow_f_plot)).transpose().tolist()
-                codes = [Path.LINETO]*len(verts)
+                verts = np.vstack(
+                    (flow_z_plot, flow_f_plot)).transpose().tolist()
+                codes = [Path.LINETO] * len(verts)
                 codes[0] = Path.MOVETO
                 path = Path(verts, codes)
                 #patch = patches.PathPatch(path, facecolor='none', lw=0.2, edgecolor='white',zorder = 2.7)
-                patch = patches.PathPatch(path, facecolor='none', lw=1.5, edgecolor='white',zorder = 1.9)
+                patch = patches.PathPatch(
+                    path, facecolor='none', lw=outline_width, edgecolor='white', zorder=1.9, alpha=0.7)
+                # patch = patches.PathPatch(
+                #     path, facecolor='none', lw=0.7, edgecolor='white', zorder=1.9, alpha=0.7)
                 ax.add_patch(patch)
-                #ax.plot(flow_z_plot, flow_f_plot, 'x',ms=2, mew=0.1, linewidth=0.5, color='k', fillstyle='none')
-
-        plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
-#                    +field_to_plot+".png")
-                    +field_to_plot+".pdf")
-        plt.draw()
+#                ax.plot(flow_z_plot, flow_f_plot, 'x', ms=2, mew=0.1,
+#                        linewidth=0.5, color='k', fillstyle='none')
 
-    #    plt.show()
-
-        plt.clf()
-        plt.close()
     finally:
-        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(np.array([x]),
-                                                                         np.array([m]))
-        print("Qabs = "+str(Qabs));
+        terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(
+            np.array([x]), np.array([m]))
+        print("Qabs = " + str(Qabs))
     #
-
-

examples/go-cc-examples.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+path=$PWD
+PROGRAM='scattnlay-example.bin'
+
+echo Compile with gcc
+rm -f $PROGRAM
+
+file=example-get-Mie.cc
+g++ -Ofast -std=c++11 $file ../src/nmie.cc ../src/nmie-applied.cc ./read-spectra.cc -lm -lrt -o $PROGRAM -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2
+
+./$PROGRAM
+# #result

src/nmie-applied.cc

@@ -75,8 +75,8 @@ namespace nmie {
         throw std::invalid_argument("Declared number of sample for Theta is not correct!");
     try {
       MultiLayerMieApplied ml_mie;
-      ml_mie.SetAllLayersSize(x);
-      ml_mie.SetAllLayersIndex(m);
+      ml_mie.SetLayersSize(x);
+      ml_mie.SetLayersIndex(m);
       ml_mie.SetAngles(Theta);
       ml_mie.SetPECLayer(pl);
       ml_mie.SetMaxTerms(nmax);
@@ -437,5 +437,27 @@ c    MM + 1  and - 1, alternately
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void MultiLayerMieApplied::RunMieCalculation() {
+    ConvertToSP();
+    MultiLayerMie::RunMieCalculation(); 
+  }
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  void MultiLayerMieApplied::GetExpanCoeffs( std::vector< std::vector<std::complex<double> > >& aln, std::vector< std::vector<std::complex<double> > >& bln, std::vector< std::vector<std::complex<double> > >& cln, std::vector< std::vector<std::complex<double> > >& dln) {
+    ConvertToSP();  // Case of call before running full Mie calculation.
+    // Calculate scattering coefficients an_ and bn_
+    calcScattCoeffs();
+    // Calculate expansion coefficients aln_,  bln_, cln_, and dln_
+    calcExpanCoeffs();
+    aln = aln_;
+    bln = bln_;
+    cln = cln_;
+    dln = dln_;
+    
+  }  // end of void MultiLayerMieApplied::GetExpanCoeffs( ...)
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
 
 }  // end of namespace nmie

src/nmie-applied.h

@@ -47,6 +47,7 @@ namespace nmie {
   class MultiLayerMieApplied : public MultiLayerMie {
     // Will throw for any error!
    public:
+    void RunMieCalculation();
     void GetFailed();
     long iformat = 0;
     bool output = true;
@@ -119,7 +120,13 @@ namespace nmie {
     std::vector<double> GetPatternEkSP();
     std::vector<double> GetPatternHkSP();
     std::vector<double> GetPatternUnpolarizedSP();
-    
+
+    void GetExpanCoeffs
+      (std::vector< std::vector<std::complex<double> > >& aln,
+       std::vector< std::vector<std::complex<double> > >& bln,
+       std::vector< std::vector<std::complex<double> > >& cln,
+       std::vector< std::vector<std::complex<double> > >& dln);
+
 
     // Output results (data file + python script to plot it with matplotlib)
     void PlotSpectra();

src/nmie.cc

@@ -1045,7 +1045,7 @@ namespace nmie {
   //**********************************************************************************//
   void MultiLayerMie::calcExpanCoeffs() {
     if (!isScaCoeffsCalc_)
-      throw std::invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
+      throw std::invalid_argument("(calcExpanCoeffs) You should calculate external coefficients first!");
 
     isExpCoeffsCalc_ = false;
 

src/nmie.h

@@ -113,7 +113,10 @@ namespace nmie {
     std::vector<double> size_param_;
     // Refractive index for all layers
     std::vector< std::complex<double> > refractive_index_;
-    // Scattering angles for scattering pattern in radians
+    // Scattering coefficients
+    std::vector<std::complex<double> > an_, bn_;
+    std::vector< std::vector<std::complex<double> > > aln_, bln_, cln_, dln_;
+    void calcExpanCoeffs();
 
   private:
     void calcNstop();
@@ -142,7 +145,6 @@ namespace nmie {
                        const double& Pi, const double& Tau, const double& n,
                        std::vector<std::complex<double> >& Mo1n, std::vector<std::complex<double> >& Me1n, 
                        std::vector<std::complex<double> >& No1n, std::vector<std::complex<double> >& Ne1n);
-    void calcExpanCoeffs();
 
     void calcField(const double Rho, const double Theta, const double Phi,
                    std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H);
@@ -151,6 +153,7 @@ namespace nmie {
     bool isScaCoeffsCalc_ = false;
     bool isMieCalculated_ = false;
 
+    // Scattering angles for scattering pattern in radians
     std::vector<double> theta_;
     // Should be -1 if there is no PEC.
     int PEC_layer_position_ = -1;
@@ -158,10 +161,7 @@ namespace nmie {
     // with calcNmax(int first_layer);
     int nmax_ = -1;
     int nmax_preset_ = -1;
-    // Scattering coefficients
-    std::vector<std::complex<double> > an_, bn_;
     std::vector< std::vector<double> > coords_;
-    std::vector< std::vector<std::complex<double> > > aln_, bln_, cln_, dln_;
     /// Store result
     double Qsca_ = 0.0, Qext_ = 0.0, Qabs_ = 0.0, Qbk_ = 0.0, Qpr_ = 0.0, asymmetry_factor_ = 0.0, albedo_ = 0.0;
     std::vector<std::vector< std::complex<double> > > E_, H_;  // {X[], Y[], Z[]}

tests/c++/speed-test-applied.cc

@@ -0,0 +1,237 @@
+//**********************************************************************************//
+//    Copyright (C) 2009-2015  Ovidio Pena <ovidio@bytesfall.com>                   //
+//    Copyright (C) 2013-2015  Konstantin Ladutenko <kostyfisik@gmail.com>          //
+//                                                                                  //
+//    This file is part of scattnlay                                                //
+//                                                                                  //
+//    This program is free software: you can redistribute it and/or modify          //
+//    it under the terms of the GNU General Public License as published by          //
+//    the Free Software Foundation, either version 3 of the License, or             //
+//    (at your option) any later version.                                           //
+//                                                                                  //
+//    This program is distributed in the hope that it will be useful,               //
+//    but WITHOUT ANY WARRANTY; without even the implied warranty of                //
+//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 //
+//    GNU General Public License for more details.                                  //
+//                                                                                  //
+//    The only additional remark is that we expect that all publications            //
+//    describing work using this software, or all commercial products               //
+//    using it, cite the following reference:                                       //
+//    [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by           //
+//        a multilayered sphere," Computer Physics Communications,                  //
+//        vol. 180, Nov. 2009, pp. 2348-2354.                                       //
+//                                                                                  //
+//    You should have received a copy of the GNU General Public License             //
+//    along with this program.  If not, see <http://www.gnu.org/licenses/>.         //
+//**********************************************************************************//
+
+#include <algorithm>
+#include <complex>
+#include <functional>
+#include <iostream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+#include <string.h>
+//sudo aptitude install libgoogle-perftools-dev
+//#include <google/heap-profiler.h>
+#include "../../src/nmie-applied.h"
+
+timespec diff(timespec start, timespec end);
+const double PI=3.14159265358979323846;
+template<class T> inline T pow2(const T value) {return value*value;}
+
+
+//***********************************************************************************//
+// This is the main function of 'scattnlay', here we read the parameters as          //
+// arguments passed to the program which should be executed with the following       //
+// syntaxis:                                                                         //
+// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment]  //
+//                                                                                   //
+// When all the parameters were correctly passed we setup the integer L (the         //
+// number of layers) and the arrays x and m, containing the size parameters and      //
+// refractive indexes of the layers, respectively and call the function nMie.        //
+// If the calculation is successful the results are printed with the following       //
+// format:                                                                           //
+//                                                                                   //
+//    * If no comment was passed:                                                    //
+//        'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                                    //
+//                                                                                   //
+//    * If a comment was passed:                                                     //
+//        'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo'                           //
+//***********************************************************************************//
+int main(int argc, char *argv[]) {
+  try {
+    std::vector<std::string> args;
+    args.assign(argv, argv + argc);
+    std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + args[0]
+			  + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] "
+			  + "[-t ti tf nt] [-c comment]\n");
+    enum mode_states {read_L, read_x, read_mr, read_mi, read_ti, read_tf, read_nt, read_comment};
+    // for (auto arg : args) std::cout<< arg <<std::endl;
+    std::string comment;
+    int has_comment = 0;
+    int i, l, L = 0;
+    std::vector<double> x, Theta;
+    std::vector<std::complex<double> > m, S1, S2;
+    double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo;
+
+    std::vector<std::complex<double> > mw, S1w, S2w;
+    double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow;
+
+    double ti = 0.0, tf = 90.0;
+    int nt = 0;    
+    if (argc < 5) throw std::invalid_argument(error_msg);
+    
+    //strcpy(comment, "");
+    // for (i = 1; i < argc; i++) {
+    int mode = -1; 
+    double tmp_mr;
+    for (auto arg : args) {
+      // For each arg in args list we detect the change of the current
+      // read mode or read the arg. The reading args algorithm works
+      // as a finite-state machine.
+
+      // Detecting new read mode (if it is a valid -key) 
+      if (arg == "-l") {
+	mode = read_L;
+	continue;
+      }
+      if (arg == "-t") {
+	if ((mode != read_x) && (mode != read_comment))
+	  throw std::invalid_argument(std::string("Unfinished layer!\n")
+							 +error_msg);
+	mode = read_ti;
+	continue;
+      }
+      if (arg == "-c") {
+	if ((mode != read_x) && (mode != read_nt))
+	  throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg);
+	mode = read_comment;
+	continue;
+      }
+      // Reading data. For invalid date the exception will be thrown
+      // with the std:: and catched in the end.
+      if (mode == read_L) {
+	L = std::stoi(arg);
+	mode = read_x;
+	continue;
+      }
+      if (mode == read_x) {
+	x.push_back(std::stod(arg));
+	mode = read_mr;
+	continue;
+      }
+      if (mode == read_mr) {
+	tmp_mr = std::stod(arg);
+	mode = read_mi;
+	continue;
+      }
+      if (mode == read_mi) {
+	m.push_back(std::complex<double>( tmp_mr,std::stod(arg) ));
+	mode = read_x;
+	continue;
+      }
+      if (mode == read_ti) {
+	ti = std::stod(arg);
+	mode = read_tf;
+	continue;
+      }
+      if (mode == read_tf) {
+	tf = std::stod(arg);
+	mode = read_nt;
+	continue;
+      }
+      if (mode == read_nt) {
+	nt = std::stoi(arg);
+        Theta.resize(nt);
+        S1.resize(nt);
+        S2.resize(nt);
+        S1w.resize(nt);
+        S2w.resize(nt);
+	continue;
+      }
+      if (mode ==  read_comment) {
+	comment = arg;
+        has_comment = 1;
+	continue;
+      }
+    }
+    if ( (x.size() != m.size()) || (L != x.size()) ) 
+      throw std::invalid_argument(std::string("Broken structure!\n")
+							 +error_msg);
+    if ( (0 == m.size()) || ( 0 == x.size()) ) 
+      throw std::invalid_argument(std::string("Empty structure!\n")
+							 +error_msg);
+    
+    if (nt < 0) {
+      printf("Error reading Theta.\n");
+      return -1;
+    } else if (nt == 1) {
+      Theta[0] = ti*PI/180.0;
+    } else {
+      for (i = 0; i < nt; i++) {
+      Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0;
+      }
+    }
+    timespec time1, time2;
+    long cpptime_nsec, best_cpp;
+    long ctime_nsec, best_c;
+    long cpptime_sec, ctime_sec;
+    long repeats = 150;
+    //HeapProfilerStart("heapprof");    
+    do {
+      clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1);
+      for (int i = 0; i<repeats; ++i) {
+    	nmie::nMieApplied(L, x, m, nt, Theta, &Qextw, &Qscaw,
+    			   &Qabsw, &Qbkw, &Qprw, &gw, &Albedow, S1w, S2w);
+      }
+      clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time2);
+      cpptime_nsec = diff(time1,time2).tv_nsec;
+      cpptime_sec = diff(time1,time2).tv_sec;
+
+      printf("-- C++ time consumed %lg sec\n", (cpptime_nsec/1e9));
+      repeats *= 10;
+    } while (cpptime_nsec < 1e8 && ctime_nsec < 1e8);
+
+        printf("\n");
+    
+    if (has_comment) {
+      printf("%6s, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", comment.c_str(), Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    } else {
+      printf("%+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e, %+.5e  \n", Qextw, Qscaw, Qabsw, Qbkw, Qprw, gw, Albedow);
+    }
+    
+    if (nt > 0) {
+      printf(" Theta,         S1.r,         S1.i,         S2.r,         S2.i\n");
+      
+      for (i = 0; i < nt; i++) {
+        printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e  \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag());
+      }
+    }
+
+  } catch( const std::invalid_argument& ia ) {
+    // Will catch if  multi_layer_mie fails or other errors.
+    std::cerr << "Invalid argument: " << ia.what() << std::endl;
+    return -1;
+  }  
+    return 0;
+}
+
+
+
+timespec diff(timespec start, timespec end)
+{
+	timespec temp;
+	if ((end.tv_nsec-start.tv_nsec)<0) {
+		temp.tv_sec = end.tv_sec-start.tv_sec-1;
+		temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec;
+	} else {
+		temp.tv_sec = end.tv_sec-start.tv_sec;
+		temp.tv_nsec = end.tv_nsec-start.tv_nsec;
+	}
+	return temp;
+}