package seed.minerva.mse; import algorithmrepository.LinearInterpolation1D; import otherSupport.Profiler; import jafama.FastMath; import oneLiners.OneLiners; import seed.minerva.ConnectionPoint; import seed.minerva.NodeImpl; import seed.minerva.StateFullNodeImpl; import seed.minerva.diagnostics.ServiceManager; import seed.minerva.neutralBeams.BeamVelocityDistributions; import seed.minerva.neutralBeams.BeamVelocityDistributionMomentSets; import seed.minerva.nodetypes.ScalarND; import seed.minerva.nodetypes.VectorND; import seed.minerva.physics.LineDefWithPerps; import seed.minerva.physics.MagneticField; import seed.minerva.physics.StokesGaussianSpectralComponents; import seed.minerva.physics.StokesSpectralComponentSetsFromMultiLine; import seed.minerva.physics.adas.AdasDataADF21_22; /** Provides lists of Stokes components of the spectrum that would be observed along a given line * from a list of velocity distributions, E and B at each point. * * Currently ignores the possible correlation between doppler shift and polarisation that would * result from the cross-component terms of the LOS parallel and LOS perp parts of the velocity distribution. * * Goes up to cubic in the stark splitting calculations (from a book). * * FIXME: Currently only works with H-alpha (not D-alpha, etc) * FIXME: The intensities were measured in pixels from a scan of the book. Eyurk! * * @author oliford */ public abstract class HAlphaMSEEmissionCore extends StateFullNodeImpl { private static final double primaryEmissionWavelength = 656.1e-9; //656.1 nm, hardcoded H-Alpha for now private static final double c = 299792458; //m s^-1 private static final int BEAM_ATOMIC_NO = 1; //Hydrogenic private static final int PLASMA_ATOMIC_NO = 1; //We assume the whole plasma is hydrogenic for now /** Component (polarisation) type definitions */ public static final int CMPT_ALL = 0; public static final int CMPT_SIGMA = 1; public static final int CMPT_PI = 2; public static final int CMPT_SIGMA_APRX_UNPOL = 3; /** Table Index definitions for starkLines 2nd dimension */ public static final int STARKINFO_POL = 0; public static final int STARKINFO_I_STAT = 1; public static final int STARKINFO_E_LINEAR = 2; public static final int STARKINFO_E_QUADRATIC = 3; public static final int STARKINFO_E_CUBIC = 4; public static final int STARKINFO_LAMBA_VARIANCE = 5; public static final int STARKINFO_NONSTAT_IDX = 6; /** Table Index definitions for non-statistical distribution corrections */ public static final int NONSTAT_NE = 0; public static final int NONSTAT_SIGMA0 = 1; public static final int NONSTAT_SIGMA1 = 2; public static final int NONSTAT_PI2 = 3; public static final int NONSTAT_PI3 = 4; public static final int NONSTAT_PI4 = 5; public static final int NONSTAT_SIGMAALL = 6; public static final int NONSTAT_PIWING = 7; /** If true, only the projection angles are really calculated * everything else is made up */ protected boolean simpleAngleCalc = true; public double neUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronDensity"); /** Wavelength splitting coefficients, up to cubic, in Hydrogen-alpha * * From figure 1.17, Chapter 17 ('The Stark Effect') from * Condon and Shortley, 'The Theory of Atomic Spectra' * * Delta (1/wavelength) = a.E + b.E^2 + c.E^3. * In the book: wavelength in cm, and E in MV cm^-1 (yes, MV, not kV like in equ10... grumble) * In the book, the equ is +/- a.E - b.E^2 +/- c.E^3 but the signs are inside the table here * * * The D-alpha ones are apparently shifted from there, but.... meh * * The intensities here were measured by hand from the graph in the book because * I spent 2 hours trying to find these numbers written down ANYWHERE, but gave up. * * 11/8/11: The 'I' in this table are the intensities when assuming a 'statistical distribution', * which is basically wrong. Alternatively nsIdx gives the index into the nonStatRatios * table for better values (which are ne dependent). * * FIXME: Probably get this from ADAS, if you ever find someone who knows which phase of the moon and what * magic incantations are required. * * * Node has no state, because it only provides a request-based method and I don't yet see any point * in caching them. * * var is effective variance due to sub-splitting of that component (used for starkLinesSimple only) * */ private static final double starkLinesFull[][] = new double[][]{ //component, I(stat), a, b, c var nsIdx { +CMPT_PI, 0.0385, +128.78, -6.715, +0.003, 0, NONSTAT_PI2, }, // (k1,k2,m) = 110 --> 010 (2(-) pi+) { +CMPT_PI, 0.1222, +193.17, -6.207, +0.088, 0, NONSTAT_PI3, }, // (k1,k2,m) = 101 --> 001 (3(0) pi+) { +CMPT_PI, 0.0892, +257.56, -6.309, +0.164, 0, NONSTAT_PI4, }, // (k1,k2,m) = 200 --> 100 (4(+) pi+) //and now all again, for the other side { -CMPT_PI, 0.0385, -128.78, -6.715, -0.003, 0, NONSTAT_PI2, }, // (k1,k2,m) = 110 --> 010 (2(-) pi-) { -CMPT_PI, 0.1222, -193.17, -6.207, -0.088, 0, NONSTAT_PI3, }, // (k1,k2,m) = 101 --> 001 (3(0) pi-) { -CMPT_PI, 0.0892, -257.56, -6.309, -0.164, 0, NONSTAT_PI4, }, // (k1,k2,m) = 200 --> 100 (4(+) pi-) //a = [0.0385, 0.1222, 0.0892]; b= [128.78, 193.17, 257.56]; //mean = 206.23352, var = 1948 // I'm not totally sure whether the sigma1 is one L circ and one R circ in different wavelen, or one of each at each wavelength. // Marteen De.Bock says he's sure that it's one of each at each wavelen, like this: { -CMPT_SIGMA, 0.0523, -64.39, -6.156, -0.085, 0, NONSTAT_SIGMA1, }, // (k1,k2,m) = 101 --> 100 (1(-) sigma+) { +CMPT_SIGMA, 0.0523, -64.39, -6.156, -0.085, 0, NONSTAT_SIGMA1, }, // (k1,k2,m) = 101 --> 100 (1(-) sigma+) { +CMPT_SIGMA, 0.0727, +0.00, -5.177, +0.000, 0, NONSTAT_SIGMA0, }, // (k1,k2,m) = 002 --> 001 (0(0) sigma+) { -CMPT_SIGMA, 0.0727, -0.00, -5.177, -0.000, 0, NONSTAT_SIGMA0, }, // (k1,k2,m) = 002 --> 001 (0(0) sigma-) { +CMPT_SIGMA, 0.0727, +0.00, -6.705, +0.000, 0, NONSTAT_SIGMA0, }, // (k1,k2,m) = 110 --> 001 (0(0) sigma+) { -CMPT_SIGMA, 0.0727, -0.00, -6.705, -0.000, 0, NONSTAT_SIGMA0, }, // (k1,k2,m) = 110 --> 001 (0(0) sigma-) { +CMPT_SIGMA, 0.0523, +64.39, -6.156, +0.085, 0, NONSTAT_SIGMA1, }, // (k1,k2,m) = 101 --> 100 (1(+) sigma+) { -CMPT_SIGMA, 0.0523, +64.39, -6.156, +0.085, 0, NONSTAT_SIGMA1, }, // (k1,k2,m) = 101 --> 100 (1(+) sigma+) //*/ //Having that many components and letting it cancel later is a bit of a waste of CPU time, so we can do this instead: //where +3 is like +/-1 but the circular stuff is incoherent instead. /* { +CMPT_SIGMA_APRX_UNPOL, 0.1046, +64.39, -6.156, +0.085, 0, , }, // (k1,k2,m) = 101 --> 100 (1(-) sigma-) { +CMPT_SIGMA_APRX_UNPOL, 0.1454, -0.00, -5.177, -0.000, 0, , }, // (k1,k2,m) = 002 --> 001 (0(0) sigma-) { +CMPT_SIGMA_APRX_UNPOL, 0.1454, -0.00, -6.705, -0.000, 0, , }, // (k1,k2,m) = 110 --> 001 (0(0) sigma-) { +CMPT_SIGMA_APRX_UNPOL, 0.1046, -64.39, -6.156, -0.085, 0, , }, // (k1,k2,m) = 101 --> 100 (1(+) sigma-) //*/ //a = [0.1046, 0.1454, 0.1454, 0.1046]; //b = [+64.39, 0, 0, -64.39]; // avg = 0, var = 1734 }; /** Primary splitting only, * used by getSimpleMSEInfo() --> calcMSESpectrum() --> getSimpleMSEInfo() */ private static final double starkLinesSimple[][] = new double[][]{ /* component I, a, b, c var nsIdx */ { +CMPT_PI, 0.25, +206.23, 0, 0, 1948, NONSTAT_PIWING, }, { +CMPT_SIGMA_APRX_UNPOL, 0.50, +0, 0, 0, 1734.7, NONSTAT_SIGMAALL, }, { -CMPT_PI, 0.25, -206.23, 0, 0, 1948, NONSTAT_PIWING, }, }; /** First/Second level splitting switch * Used by getStokesComponents() --> calcMSESpectrum() --> fillStokesComponents() */ //private static double starkLines[][] = starkLinesSimple; private static double starkLines[][] = starkLinesFull; /** Fractional component intensities for non-statistical distribution * Taken from O Marchuk et al. "Collisional excitation and emission of Hα * Stark multiplet in fusion plasmas" 2010 J. Phys.B 43 011002 * http://dx.doi.org/10.1088/0953-4075/43/1/011002 * and his accompanying talk at IPP Greifswald 05/05/2011. * * I think these are valid for Te~3keV, B=1T, E ~perp to v * * First entry is log10(ne/m^-3) and MUST be regular */ public final static double nonStatRatios[][] = new double[][]{ { 17.000, 17.300, 17.600, 17.900, 18.200, 18.500, 18.800, 19.100, 19.400, 19.700, 20.000, 20.300, 20.600, 20.900 }, //ne { 0.045, 0.046, 0.048, 0.052, 0.057, 0.062, 0.069, 0.074, 0.079, 0.083, 0.086, 0.088, 0.091, 0.094 }, //sigma0 { 0.075, 0.075, 0.074, 0.074, 0.073, 0.071, 0.070, 0.070, 0.069, 0.068, 0.067, 0.067, 0.066, 0.065 }, //sigma1 { 0.050, 0.049, 0.047, 0.044, 0.043, 0.041, 0.039, 0.038, 0.037, 0.037, 0.038, 0.038, 0.038, 0.038 }, //pi2 { 0.121, 0.121, 0.121, 0.121, 0.121, 0.122, 0.122, 0.122, 0.122, 0.122, 0.122, 0.122, 0.122, 0.122 }, //pi3 { 0.164, 0.163, 0.161, 0.157, 0.150, 0.142, 0.132, 0.122, 0.114, 0.106, 0.102, 0.096, 0.091, 0.086 }, //pi4 { 0.331, 0.336, 0.342, 0.355, 0.372, 0.392, 0.415, 0.436, 0.453, 0.466, 0.476, 0.487, 0.496, 0.507 }, //sigmaAll { 0.335, 0.333, 0.329, 0.323, 0.314, 0.304, 0.293, 0.282, 0.274, 0.267, 0.262, 0.257, 0.252, 0.247 }, //piWing }; /** Parent from beams that gives a list of multivariate Gaussian beam particle velocity * distributions for each requested point */ BeamVelocityDistributionMomentSets neutralBeamVelocities; /** Magnetic field, B_{XYZ}(x,y,z) in global/lab rest frame */ MagneticField magneticField; /** Electric field, E_{XYZ}(x,y,z) in global/lab rest frame */ VectorND electricField; /** Electron density */ ScalarND electronDensity; /** Hack to arbitarily rotate the B field a bit. * If set to ±inf, only toroidal field component is kept nonzero*/ public double fieldHackeryRotation = 0; /** Force output of radiation in a single component type: one of CMPT_TYPE_* */ public int forceSingleComponent = CMPT_ALL; private boolean nonStatIntensities = false; public HAlphaMSEEmissionCore() { addConnectionPoint(new ConnectionPoint("neutralBeamVelocities", BeamVelocityDistributionMomentSets.class, false, getField("neutralBeamVelocities"))); addConnectionPoint(new ConnectionPoint("magneticField", MagneticField.class, false, getField("magneticField"))); addConnectionPoint(new ConnectionPoint("electricField", VectorND.class, true, getField("electricField"))); addConnectionPoint(new ConnectionPoint("electronDensity", ScalarND.class, true, getField("electronDensity"))); } public static double B2[] = new double[3]; /** Guts of the MSE calculation * Expects all vectors to be in LUR coordinates (line of sight, 'up', 'right'): * E, B, v (velDists.mean, .covar) * * @param posIdx List of indices into velDists * @param k0 starting index along line to do * @param k1 ending index along line to do * */ protected final Object calcMSESpectrum(BeamVelocityDistributions velDists, boolean fullStokes, int posIdx[], double losE[][], double losB[][], double ne[], int k0, int k1) { final int nLines = starkLines.length; //we need to be able to evaluate the effective beam emission coefficients AdasDataADF21_22 qee = ServiceManager.getInstance().getADASReader().getADF21_22DataSet(true, BEAM_ATOMIC_NO, 0, //beamAtomNumber, beamAtomCharge, PLASMA_ATOMIC_NO, 1); //plasmaAtomNumber, plasmaAtomCharge); StokesGaussianSpectralComponents stokes = null; SimpleMSEPolarisationInfos mseInfos = null; // for each velcoity component, we have each stark component if(fullStokes) stokes = new StokesGaussianSpectralComponents(velDists.getNSets(), velDists.maxDistribs*nLines); else //or just 1 set of simple info for each velocity component mseInfos = new SimpleMSEPolarisationInfos(velDists.getNSets(), velDists.maxDistribs); for(int k=k0; k <= k1; k++){ //for each position along the LOS we are doing (offset by k0)... int iP = posIdx[k]; //index onto line if(Double.isNaN(losB[0][k]) || ne[k] == 0) //no field data out here continue; int nDists = velDists.getN(iP); if(fullStokes) stokes.nStokes[iP] = nDists*nLines; else //or just 1 set of simple info for each velocity component mseInfos.nComponents[iP] = nDists; for(int iD=0; iD < nDists; iD++){ //for each gaussian velocity distribution component int iPiD = iP*velDists.maxDistribs+iD; if(velDists.atomicNumber[iPiD] != BEAM_ATOMIC_NO) throw new IllegalArgumentException("HAlphaMSEEmission: found non-hydrogren (atomic no = " + velDists.atomicNumber[iPiD] + ") neutral beam component."); // The ADAS effective emission coefficient qee, is in m^3 s^-1 // The reaction cross-section, in m^2, is qe = qee / beam velocity // The reaction rate (and photon generation rate), in m^-3 s-^1, is R = qe * ne * neutral flux // neutral flux = neutral density * beam velocity // so we end up cancelling the beam velocity to get R = qe * ne * neutral density double photonRate = qee.evalCoeff(velDists.meanEnergy[iPiD], ne[k]*neUnit) * (ne[k]*neUnit) * velDists.neutralDensity[iPiD]; //double photonRate = velDists[i][j].neutralDensity; double meanLUR[] = velDists.getMean(iP, iD); double covarLUR[] = velDists.getCovar(iP, iD); // calc doppler shift from velocity parallel to los (observer assumed sitting at the line.P0()) double lambda0 = ((1 + meanLUR[0]/c) * primaryEmissionWavelength); //there is also a variance in wavelength due to the variance in beam parallel velocity // FIXME: we are ignoring the correlation here, because it's too hard to think about double varLambda = covarLUR[0] * primaryEmissionWavelength*primaryEmissionWavelength / (c*c); //For nonstatistical distributions, we have to interpolate the table with ne, some of which can be done now int nonStatIntensityIdx=-1; double nonStatIntensityFrac=0; if(nonStatIntensities){ double d = FastMath.log10(ne[k]*neUnit); nonStatIntensityFrac = ((d - nonStatRatios[0][0])*(nonStatRatios[0].length-1)/(nonStatRatios[0][nonStatRatios[0].length-1]-nonStatRatios[0][0])); nonStatIntensityIdx = (int)nonStatIntensityFrac; nonStatIntensityFrac -= nonStatIntensityIdx; if(nonStatIntensityIdx < 0){ nonStatIntensityIdx = 0; nonStatIntensityFrac = 0; } if(nonStatIntensityIdx >= nonStatRatios[0].length){ nonStatIntensityIdx = nonStatRatios[0].length-1; nonStatIntensityFrac = 1; } } //get the (average) local E field in the beam particles rest frame E' = E + v ^ B double localLosE[] = new double[] { +meanLUR[1]*losB[2][k] - meanLUR[2]*losB[1][k], -meanLUR[0]*losB[2][k] + meanLUR[2]*losB[0][k], +meanLUR[0]*losB[1][k] - meanLUR[1]*losB[0][k] }; if(losE != null){ //FIXME: Check if we can really just add this like this. localLosE[0] += losE[0][k]; localLosE[1] += losE[1][k]; localLosE[2] += losE[2][k]; } //we need the magnitude of the E field double magE = Math.sqrt(localLosE[0]*localLosE[0] + localLosE[1]*localLosE[1] + localLosE[2]*localLosE[2]); double magE_MVperCM = magE / 1e6 / 100; //(in MV/cm, otherwise the coefficients are silly orders) //we can do the polarisation trig out here double cosGamma = localLosE[0] / magE; // angle between E and LOS (gamma, in Howard2008 notation) double cosGammaSq = cosGamma*cosGamma; // trig avoidance double sinGammaSq = 1 - cosGammaSq; // beta = angle of electric field as angle from LOS 'up' definition // in the +ve R and +ve U direction (i.e. clockwise looking along parallel to LOS) // PI is polarised parallel to E, and sigma at +ve 90' further around in the same // direction (as beta) //tan(beta) = localLosE[2] / localLosE[1], so.... (trig avoidance) double EruLenSq = (localLosE[1]* localLosE[1] + localLosE[2]* localLosE[2]); double cos2Beta = (localLosE[1]* localLosE[1] - localLosE[2]* localLosE[2]) / EruLenSq; double sin2Beta = 2 * localLosE[1] * localLosE[2] / EruLenSq; if(EruLenSq <= 0){ cos2Beta = 1; sin2Beta = 0; } //entirely unpolarised, but we need to be able todo maths if(fullStokes) fillStokesComponents(stokes, iP, iD, magE_MVperCM, cosGamma, cosGammaSq, sinGammaSq, cos2Beta, sin2Beta, lambda0, varLambda, photonRate, nonStatIntensityIdx, nonStatIntensityFrac); else fillSimpleMSEInfo(mseInfos, iP, iD, magE_MVperCM, cosGamma, cosGammaSq, sinGammaSq, cos2Beta, sin2Beta, lambda0, varLambda, photonRate, nonStatIntensityIdx, nonStatIntensityFrac, magE); } } return fullStokes ? stokes : mseInfos; } /** Internal for calcMSESpectrum() when called by getSimpleMSEInfo() */ private final void fillSimpleMSEInfo(SimpleMSEPolarisationInfos mseInfo, int iP, int iC, double magE_MVperCM, double cosGamma, double cosGammaSq, double sinGammaSq, double cos2Beta, double sin2Beta, double lambda0, double varLambda, double photonRate, int nsiIdx, double nsiFrac, double magE){ int iPiC = iP * mseInfo.maxComponents + iC; //calc average shift of PI commponents (linear only) double deltaPerCMPI = starkLinesSimple[0][STARKINFO_E_LINEAR] * magE_MVperCM; // Calc shifted wavelength mseInfo.lambda0[iPiC] = lambda0; mseInfo.deltaLambdaPI[iPiC] = lambda0 * (1.0 - 1.0 / (1.0 + 100 * lambda0 * deltaPerCMPI)); // Divide intensity amongst the stark components //FIXME: This needs the non-statistical stuff adding mseInfo.Ipi[iPiC] = photonRate * sinGammaSq * starkLinesSimple[0][STARKINFO_I_STAT]; mseInfo.IsigmaPol[iPiC] = photonRate * sinGammaSq * starkLinesSimple[1][STARKINFO_I_STAT]; mseInfo.IsigmaUnpol[iPiC] = photonRate * (2 * cosGammaSq) * starkLinesSimple[1][STARKINFO_I_STAT]; //and variance in that due to 'variance' of sub-level splitting double varLambdaSubLevelPi = lambda0*lambda0*lambda0*lambda0 * starkLinesSimple[0][STARKINFO_LAMBA_VARIANCE] * 10000 * magE_MVperCM * magE_MVperCM; double varLambdaSubLevelSigma = lambda0*lambda0*lambda0*lambda0 * starkLinesSimple[1][STARKINFO_LAMBA_VARIANCE] * 10000 * magE_MVperCM * magE_MVperCM; // Variance in wavelength due to variance in doppler shift and variance due to sub-level splitting // This ignores the reshaping of that distribution due to the stark shift, but it will be VERY small. mseInfo.lambdaVarSigma[iPiC] = varLambda + varLambdaSubLevelSigma; mseInfo.lambdaVarPi[iPiC] = varLambda + varLambdaSubLevelPi; //mseInfo.anglePiToR[iPiC] = FastMath.atan2(sin2Beta, cos2Beta) / 2; mseInfo.sin2PiFromU[iPiC] = sin2Beta; //trig avoidance mseInfo.cos2PiFromU[iPiC] = cos2Beta; //trig avoidance mseInfo.magE[iPiC] = magE; } /** Internal for calcMSESpectrum() when called by getStokesComponents() */ private final void fillStokesComponents(StokesGaussianSpectralComponents stokesVecs, int iP, int j, double magE_MVperCM, double cosGamma, double cosGammaSq, double sinGammaSq, double cos2Beta, double sin2Beta, double lambda0, double varLambda, double photonRate, int nsiIdx, double nsiFrac){ final int nLines = starkLines.length; for(int iSC = 0; iSC < nLines; iSC++){ //for each stark splitting component int iS = j*nLines + iSC; int iPiS = iP * stokesVecs.maxStokes + iS; //override for only outputting single component types if(forceSingleComponent != CMPT_ALL && forceSingleComponent != Math.abs(starkLines[iSC][STARKINFO_POL])) continue; double deltaPerCM = starkLines[iSC][STARKINFO_E_LINEAR] * magE_MVperCM + starkLines[iSC][STARKINFO_E_QUADRATIC] * magE_MVperCM * magE_MVperCM + starkLines[iSC][STARKINFO_E_CUBIC] * magE_MVperCM * magE_MVperCM * magE_MVperCM; // Calc shifted wavelength stokesVecs.lambda0[iPiS] = 1.0 / ((1.0 / lambda0) + (deltaPerCM*100)); //and variance in that due to 'variance' of sub-level splitting double varLambdaSubLevel = stokesVecs.lambda0[iPiS]*stokesVecs.lambda0[iPiS] * stokesVecs.lambda0[iPiS]*stokesVecs.lambda0[iPiS] * starkLines[iSC][STARKINFO_LAMBA_VARIANCE] * 10000 * magE_MVperCM * magE_MVperCM; if(nonStatIntensities){ int k = (int)starkLines[iSC][STARKINFO_NONSTAT_IDX]; //interpolate the table stokesVecs.I[iPiS] = photonRate * ((1-nsiFrac) * nonStatRatios[k][nsiIdx] + nsiFrac * nonStatRatios[k][nsiIdx+1]); }else{ // Divide intensity amongst the stark components stokesVecs.I[iPiS] = photonRate * starkLines[iSC][STARKINFO_I_STAT]; } // Variance in wavelength due to variance in doppler shift and variance due to sub-level splitting // This ignores the reshaping of that distribution due to the stark shift, but it will be VERY small. stokesVecs.lambdaVar[iPiS] = varLambda + varLambdaSubLevel; //and now the polarisation /* The stokes vectors are different for the sigma and PI components * and also for the +/- sigma components. * * For Starke, the +/- difference cancels but AFAIK this is then the same for * also Zeeman splitting so we can generalise this node/class later on :) * * In fact the stokes formulas are FOR Zeeman splitting and come from Howard2008 * but I'm not sure what he's trying to say is and isn't valid for MSE. * * The result from this 'appears' to make sense, but more testing is needed. * * Ho hum. * * FIXME: No attempt so far to even represent the distribution in Stokes vecs, * due to distribution in v. Actually isn't just a reduction in the * degree of polarisation???? - needs thinking through carefully * * The full stokes vectors here are: * sigma = I0 * [(1 + cosGammaSq), -sinGammaSq cos2Beta, -sinGammaSq sin2Beta, 0] * sigma = I0 * (1 + cosGammaSq) * [ 1, * -sinGammaSq cos2Beta / (1 + cosGammaSq), * -sinGammaSq sin2Beta / (1 + cosGammaSq), * 0] * * pi = I0 * sinGammaSq * [ 1, cos2Beta, sin2Beta, 0 ] * * For gamma = 0, the intensity for sigma becomes 2*I0, but pi becomes 0. * Since the I0 already includes the component intensity coefficients from the table * and sum of all the sigma components = 0.5, and sum of all the pi components = 0.5 * this will give, for gamma=0, Isig + Ipi = 1 * * but we want to write in I * [s0/I, s1/I, s2/I ] * */ if(starkLines[iSC][STARKINFO_POL] == -CMPT_SIGMA || starkLines[iSC][STARKINFO_POL] == +CMPT_SIGMA || starkLines[iSC][STARKINFO_POL] == +CMPT_SIGMA_APRX_UNPOL) { //sigma + and - summed (CPU time hack, see table) double Ifact = (1 + cosGammaSq); stokesVecs.I[iPiS] *= Ifact; stokesVecs.stokes[iPiS*3+0] = - sinGammaSq * cos2Beta / Ifact; stokesVecs.stokes[iPiS*3+1] = - sinGammaSq * sin2Beta / Ifact; stokesVecs.stokes[iPiS*3+2] = (starkLines[iSC][STARKINFO_POL] == CMPT_SIGMA_APRX_UNPOL) ? 0 : (-starkLines[iSC][STARKINFO_POL] * 2 * cosGamma / Ifact); //for stark, the +/- on this last one will cancel if all the components get summed later } else { // +/-2 = pi+/- stokesVecs.I[iPiS] *= sinGammaSq; stokesVecs.stokes[iPiS*3+0] = //sinGammaSq * cos2Beta / Ifact; cos2Beta; stokesVecs.stokes[iPiS*3+1] = //sinGammaSq * sin2Beta / Ifact; sin2Beta; stokesVecs.stokes[iPiS*3+2] = 0; } stokesVecs.tagID[iPiS] = (int)starkLines[iSC][STARKINFO_POL]; } } /** Choose between simple and full splitting modes */ public void setSimpleSplittingMode(boolean enableSimpleSplitting){ if(enableSimpleSplitting){ starkLines = starkLinesSimple; }else{ starkLines = starkLinesFull; } } /** Force outputting of only only component type: one of CMPT_* */ public void setForceSingleComponent(int forceSingleComponent){ this.forceSingleComponent = forceSingleComponent; } public int getForceSingleComponent(){ return this.forceSingleComponent; } }