// Header file containing the definition of the 'mie' routines. These // routines compute tables of scattering and extinction // cross-sections, albedo functions, and entries of the scattering // matrix for an average dust grain. This average dust grain is given // by integration over the defined dust composition and grain size // distribution. The code for these calculations is taken from Bohren // & Huffman (1983). It is based on solving the Maxwell equations for // an electromagnetic wave scattered off a spherical grain. We // implemented three tables of optical data for graphite and silicate // grains. The non-isotropic optical properties of graphite crystals // are considered by using two sorts of grains, 'GraphOrt' and // 'GraphPar', for electromagnetic waves traveling perpendicular and // parallel to the crystal's symmetry axis. The dielectric functions // were taken from Draine & Lee (1984). The computed Mie data is // stored in a look-up table of the object 'mie'. // // References: // Bohren, C. F. & Huffman, D. R. 1983, New York: Wiley // Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89 // // STOKES, version 1.0, Nov 2004 #ifndef mie_h #define mie_h #include "specmath-stokes-v1.0.h" // definition of the angular resolution in theta for the Mie data const int NANG = 91; // maximum iteration depth of the D-function const int DMax = 3000; // maximum number of entries in the dielectric function tables const int DataMax = 71; // maximum number of dust components for the Mie data table const int SubNum = 3; // maximum number of entries in the theta sampling list const int ThetaAcc = 1800; // defines an entry of the dielectric function tables struct point { double lam; double eps; }; // real part of the dielectric function for 'GraphPar' const point GraphParRe[DataMax] = {{1000.40, -0.90044}, {1015.12, -0.90044}, {1030.28, -0.83588}, {1043.93, -0.60992}, {1062.01, 0.55216}, {1082.85, 2.263}, {1104.52, 4.13524}, {1122.50, 4.8454}, {1145.81, 5.55556}, {1162.71, 5.9752}, {1182.64, 6.20116}, {1200.66, 6.1366}, {1232.85, 5.74924}, {1299.63, 4.90996}, {1370.65, 4.32892}, {1453.70, 3.84472}, {1543.13, 3.42508}, {1654.15, 3.03772}, {1765.28, 2.9086}, {1866.81, 2.7472}, {1952.68, 2.55352}, {2062.14, 2.35984}, {2117.54, 2.263}, {2202.04, 2.13388}, {2312.79, 2.03704}, {2393.03, 1.97248}, {2479.04, 1.97248}, {2524.40, 2.00476}, {2595.65, 2.03704}, {2671.03, 2.13388}, {2778.63, 2.32756}, {2895.26, 2.5858}, {2989.37, 2.97316}, {3089.80, 3.29596}, {3178.80, 3.5542}, {3292.61, 3.78016}, {3414.86, 4.0384}, {3546.55, 4.16752}, {3713.63, 4.3612}, {3842.94, 4.58716}, {3981.58, 4.71628}, {4099.91, 4.87768}, {4324.85, 5.10364}, {4501.24, 5.26504}, {4653.06, 5.26504}, {4815.49, 5.29732}, {4989.66, 5.29732}, {5081.56, 5.3296}, {5225.94, 5.36188}, {5378.75, 5.36188}, {5485.70, 5.39416}, {5596.98, 5.39416}, {5772.63, 5.36188}, {5895.99, 5.36188}, {6091.25, 5.3296}, {6228.76, 5.26504}, {6372.63, 5.23276}, {6601.33, 5.23276}, {6763.15, 5.20048}, {6933.10, 5.20048}, {7204.66, 5.1682}, {7397.84, 5.1682}, {7707.85, 5.13592}, {7929.37, 5.10364}, {8163.99, 5.07136}, {8412.93, 5.07136}, {8543.18, 5.03908}, {8816.17, 5.0068}, {9107.17, 5.0068}, {9418.05, 5.0068}, {9926.30, 5.0068}}; // imaginary part of the dielectric function for 'GraphPar' const point GraphParIm[DataMax] = {{993.61, 4.85085}, {1030.79, 5.78115}, {1052.45, 6.44565}, {1079.25, 8.37270}, {1105.22, 8.93753}, {1125.54, 7.84110}, {1173.49, 4.38570}, {1233.93, 2.39220}, {1304.02, 1.46190}, {1372.22, 1.12965}, {1451.77, 0.83063}, {1545.42, 0.73095}, {1647.06, 0.83063}, {1768.66, 0.86385}, {1922.98, 0.86385}, {2075.07, 0.99675}, {2147.87, 1.16288}, {2329.49, 1.52835}, {2476.05, 1.99350}, {2532.97, 2.25930}, {2707.21, 2.72445}, {2832.46, 2.99025}, {2938.19, 3.15638}, {3086.32, 3.25605}, {3269.45, 3.18960}, {3475.70, 2.95703}, {3660.42, 2.85735}, {3838.95, 2.82413}, {3977.52, 2.69123}, {4095.79, 2.59155}, {4157.60, 2.52510}, {4354.76, 2.32575}, {4460.52, 2.09318}, {4728.48, 1.89383}, {4985.18, 1.62803}, {5221.39, 1.42868}, {5481.10, 1.19610}, {5708.24, 1.09643}, {5955.02, 1.02998}, {6224.10, 0.93030}, {6367.97, 0.86385}, {6596.70, 0.79740}, {6758.53, 0.73095}, {7016.74, 0.69773}, {7200.13, 0.66450}, {7393.36, 0.63128}, {7597.25, 0.59805}, {7812.70, 0.56483}, {8040.73, 0.53160}, {8282.47, 0.49838}, {8539.20, 0.43193}, {8812.35, 0.39870}, {9103.55, 0.36548}, {9414.66, 0.33225}, {9747.78, 0.29903}, {10105.34, 0.26580}}; // real part of the dielectric function for 'GraphOrt' const point GraphOrtRe[DataMax] = {{1000.40, 2.84404}, {1015.12, 3.26368}, {1030.28, 3.74788}, {1043.93, 4.13524}, {1062.01, 4.29664}, {1082.85, 4.32892}, {1104.52, 4.16752}, {1122.50, 4.0384}, {1145.81, 3.84472}, {1162.71, 3.68332}, {1182.64, 3.45736}, {1200.66, 3.26368}, {1232.85, 3.10228}, {1299.63, 2.65036}, {1370.65, 2.13388}, {1453.70, 1.68196}, {1543.13, 1.26232}, {1654.15, 0.77812}, {1765.28, 0.35848}, {1866.81, -0.12572}, {1952.68, -0.54536}, {2062.14, -1.09412}, {2117.54, -1.51376}, {2202.04, -2.12708}, {2312.79, -2.86952}, {2393.03, -3.386}, {2479.04, -3.67652}, {2524.40, -3.57968}, {2595.65, -3.25688}, {2671.03, -2.64356}, {2778.63, -1.25552}, {2895.26, 0.10024}, {2989.37, 1.26232}, {3089.80, 2.16616}, {3178.80, 2.71492}, {3292.61, 3.13456}, {3414.86, 3.5542}, {3546.55, 3.74788}, {3713.63, 3.97384}, {3842.94, 4.16752}, {3981.58, 4.29664}, {4099.91, 4.39348}, {4324.85, 4.58716}, {4501.24, 4.684}, {4653.06, 4.74856}, {4815.49, 4.81312}, {4989.66, 4.90996}, {5081.56, 4.97452}, {5225.94, 5.0068}, {5378.75, 5.03908}, {5485.70, 5.07136}, {5596.98, 5.10364}, {5772.63, 5.13592}, {5895.99, 5.1682}, {6091.25, 5.23276}, {6228.76, 5.26504}, {6372.63, 5.29732}, {6601.33, 5.3296}, {6763.15, 5.36188}, {6933.10, 5.39416}, {7204.66, 5.45872}, {7397.84, 5.55556}, {7707.85, 5.6524}, {7929.37, 5.71696}, {8163.99, 5.78152}, {8412.93, 5.87836}, {8543.18, 5.94292}, {8816.17, 6.03976}, {9107.17, 6.16888}, {9418.05, 6.298}, {9926.30, 6.55624}}; // imaginary part of the dielectric function for 'GraphOrt' const point GraphOrtIm[DataMax] = {{993.61, 5.44890}, {1030.79, 4.58505}, {1052.45, 4.11990}, {1079.25, 2.99025}, {1105.22, 2.29253}, {1125.54, 1.92705}, {1173.49, 1.42868}, {1233.93, 1.02998}, {1304.02, 0.76418}, {1372.22, 0.66450}, {1451.77, 0.66450}, {1545.42, 0.73095}, {1647.06, 0.86385}, {1768.66, 1.06320}, {1922.98, 1.36223}, {2075.07, 1.76093}, {2147.87, 2.19285}, {2329.49, 3.52185}, {2476.05, 5.58180}, {2532.97, 6.97725}, {2707.21, 9.36945}, {2832.46, 10.16685}, {2938.19, 10.20008}, {3086.32, 10.00073}, {3269.45, 9.50235}, {3475.70, 9.03720}, {3660.42, 8.77140}, {3838.95, 8.63850}, {3977.52, 8.60528}, {4095.79, 8.57205}, {4157.60, 8.57205}, {4354.76, 8.43915}, {4460.52, 8.50560}, {4728.48, 8.60528}, {4985.18, 8.67173}, {5221.39, 8.70495}, {5481.10, 8.77140}, {5708.24, 8.97075}, {5955.02, 9.17010}, {6224.10, 9.40268}, {6367.97, 9.50235}, {6596.70, 9.63525}, {6758.53, 9.73493}, {7016.74, 9.83460}, {7200.13, 10.03395}, {7393.36, 10.29975}, {7597.25, 10.53233}, {7812.70, 10.76490}, {8040.73, 10.96425}, {8282.47, 11.16360}, {8539.20, 11.46263}, {8812.35, 11.66198}, {9103.55, 11.99423}, {9414.66, 12.55905}, {9747.78, 12.72518}, {10105.34, 12.89130}}; // real part of the dielectric function for 'Astronomical Silicate' const point SilicateRe[DataMax] = {{996.62, 1.823104}, {1067.33, 2.448945}, {1143.06, 2.829892}, {1224.16, 3.156418}, {1311.01, 3.319681}, {1380.17, 3.292471}, {1452.98, 3.700628}, {1556.07, 5.006732}, {1666.47, 4.598575}, {1784.71, 4.244838}, {1911.33, 3.972733}, {2082.32, 3.755049}, {2230.06, 3.537365}, {2471.56, 3.346892}, {2739.20, 3.238050}, {3035.83, 3.156418}, {3665.58, 3.074787}, {4204.18, 3.047576}, {4821.92, 3.020366}, {5344.09, 3.020366}, {6129.31, 3.020366}, {7029.91, 2.993155}, {7528.68, 2.993155}, {8343.97, 2.993155}, {9090.42, 2.993155}, {9903.66, 2.993155}, {10606.32, 2.993155}, {11555.16, 2.965945}, {12375.00, 2.965945} }; // imaginary part of the dielectric function for 'Astronomical Silicate' const point SilicateIm[DataMax] = {{979.23, 2.79382}, {1013.61, 3.04050}, {1049.21, 3.22992}, {1104.95, 2.93222}, {1163.66, 2.89699}, {1246.82, 2.59838}, {1313.07, 2.47574}, {1382.83, 2.41661}, {1481.65, 3.30895}, {1507.44, 3.00397}, {1560.38, 1.68161}, {1643.28, 0.85459}, {1760.72, 0.53986}, {1886.54, 0.30961}, {1986.78, 0.22611}, {2056.55, 0.16514}, {2241.87, 0.10183}, {2360.98, 0.09470}, {2486.43, 0.09940}, {3058.49, 0.09940}, {3697.81, 0.09940}, {4101.19, 0.10060}, {4708.31, 0.10060}, {5405.29, 0.10060}, {6205.45, 0.10060}, {7124.06, 0.10183}, {8038.74, 0.10183}, {9228.74, 0.10307}, {11352.04, 0.10432}}; // definition of the object for the Mie data table class Mie { public: // initialize the Mie data table by computing cross-sections, // albedos, and entries of the scattering matrix void InitMie(double& RadMin,double& RadMax,int& RadNum, int& LamNum,double Comp[SubNum],double& RadAlpha, char dustfile[80]); // return the total extinction cross-section of the average grain at // wavelength lambda double GetCext(double lambda); // return the albedo of the average grain at wavelength lambda double GetAlbedo(double lambda); // return the scattering angle corresponding to a random number r // and wavelength lambda double GetTheta(double lambda,double r); // return entry S11 of the Mueller matrix at wavelength lambda and // scattering angle theta double GetS11(double lambda,double theta); // return entry S12 of the Mueller matrix at wavelength lambda and // scattering angle theta double GetS12(double lambda,double theta); // return entry S33 of the Mueller matrix at wavelength lambda and // scattering angle theta double GetS33(double lambda,double theta); // return entry S34 of the Mueller matrix at wavelength lambda and // scattering angle theta double GetS34(double lambda,double theta); // set S11, S12, S33, and S34 to the entries of the Mueller matrix // at the wavelength lambda and scattering angle theta void GetMueller(double lambda,double theta,double& S11,double& S12, double& S33,double& S34); private: // define an entry in the Mie data table struct Entry { double Cext; double Albedo; double Theta[ThetaAcc]; double S11[2*NANG-1],S12[2*NANG-1],S33[2*NANG-1],S34[2*NANG-1]; }; // define type of a pointer on an entry of the Mie data table typedef Entry* EntryPtr; // declare a pointer to entries of the Mie data table, the dynamical // array over values of lambda will be defined later in the program EntryPtr MD; // declare minimum and maximum wavelength and width of one step // in lambda double GrainLamMin,GrainLamMax,GrainLamDiff; // declare number of wavelengths considered int GrainLamNum; }; #endif