package seed.minerva.polarimetry; import oneLiners.OneLiners; import oneLiners.AsciiMatrixFile; import seed.minerva.ConnectionPoint; import seed.minerva.MinervaSettings; import seed.minerva.NodeImpl; import seed.minerva.StateFullNodeImpl; import seed.minerva.diagnostics.Diagnostic; import seed.minerva.diagnostics.LineOfSightSystem; import seed.minerva.diagnostics.ServiceManager; import seed.minerva.nodetypes.*; import seed.minerva.physics.ElectronDensity; import seed.minerva.physics.MagneticField; import sun.reflect.generics.reflectiveObjects.NotImplementedException; import algorithmrepository.Algorithms; /**
* Based on: * "Review of plasma polarimetry" S. E. Segre * * Stokes representation: * s_0 = cos 2chi . cos 2psi * s_1 = cos 2chi . sin 2psi * s_2 = sin 2chi * * Final stokes from initial: * sF = M . sI * * Equations 6.15, 6.6, 6.5, 6.8: * M = / 0 -W3 W2 \ * | W3 0 -W1 | * \ -W2 W1 0 / * * W1 = c1l3 int[ ne . (Bx^2 - By^2) dz ] * W2 = c1l3 int[ ne . Bx . By dz ] * W3 = c3l2 int[ ne . Bz dz ] * * The 'analytical' prediction is based on these: * * sF_0 = - W3 sI_1 + W2 sI_2 * sF_1 = W3 sI_0 - W1 sI_2 * sF_2 = -W2 sI_0 + W1 sI_1 * * It assumes that chi0 = 0: * (1) cos 2psiF = - W3 . sin 2psi0 * (2) sin 2psiF = W3 . cos 2psi0 * (3) sin 2chiF = -W2 . cos 2psi0 + W1 sin 2psi0 * * For psi0 ~ 0: * DPsi = 1/2 W3 = 1/2 c3l2 int[ ne . Bz dz ] * 2 chiF = - 1/2 W2 = - 1/2 c1l3 int[ ne . Bx . By dz ] * * For psi0 ~ 45': * cos (90' + 2Dpsi) = -W3, so: * DPsi = - 1/2 W3 = - 1/2 c3l2 int[ ne . Bz dz ] * 2 chiF = 1/2 W1 = 1/2 c1l3 int[ ne . (Bx^2 - By^2) dz ] * (same as before, but -ve?) * * The coordinates assumed here have z // with los, y perp to both los and toroidal dir and x ^ y = z. * So for any LOS in the poloidal plane, Bx = Bphi and By ~ 0, so chiF = 0 if psi0 = 0. * * Final analytic is therefor: * DPsi = 1/2 c3l2 int[ ne . Bz dz ] * chiF = 0 for psi0 ~ 0 * chiF = 1/2 c1l3 int[ ne . Btor^2 dz ] for psi0 ~ 45' * * Wasn't that fun. * **/ public class PolarimeterPredictionNode extends StateFullNodeImpl { private static final double d2r = Math.PI / 180; private static final double r2d = 180 / Math.PI; private static final double c = 299792458; private static final double e = 1.60217646e-19; //electron charge, not exponential! private static final double me = 9.10938188e-31; private static final double ep0 = 8.85418782e-12; private static final double c1 = e*e*e*e / (16 * Math.PI*Math.PI*Math.PI * ep0 * me*me*me * c*c*c*c); private static final double c3 = e*e*e / ( 4 * Math.PI*Math.PI * ep0 * me*me * c*c*c); public final static int PLASMA_MODEL_COLD = 1; public final static int PLASMA_MODEL_WARM_SEGRE = 2; public final static int PLASMA_MODEL_WARM_MIRNOV = 3; public final static int PLASMA_MODEL_WARM = PLASMA_MODEL_WARM_MIRNOV; // Parents ScalarND electronDensity, electronTemperature; LineOfSightSystem los; IntegerArray dPsiStatus, chiStatus; MagneticField magneticField; /** Angle between initial linear polarisation's E vector and the direction perp to LOS which is * most parallel to +ve toroidal dir */ DoubleArray psi0, chi0; DoubleValue laserWavelength; //laser angular freq BooleanArray returnPass; // Properties private int predEnableRequest[]; /* Diagnotic.ON, Diagnostic.OFF or Diagnostic.AUTO to take from channelOK */ private int plasmaModel = PLASMA_MODEL_COLD; // State private int nChannels; private double dPsi[]; private double chi[]; private boolean doublePass[]; //laser angular freq double w = Double.NaN; /** whether or not we should calc DPsi and chi - dependant on prediction enable and psiStatus/chiStatus */ private boolean calcDPsi[], calcChi[]; /** Perform the analytical integrals instead of the full evolution evaluation. * The analytical version expects the initial polarisation to have chi0 ~ 0 * and psi0 ~ 0 or psi0 ~ 45 degrees. */ private boolean approxAnalytic = false; private boolean psiIs45Degrees[]; /** points along LOS [channel][{x,y,z}][point] */ private double losPoints[][][]; /** dists along LOS */ private double[][] losDists; /** unit vector along LOS [channel][{x,y,z}] */ private double uPara[][]; /** unit vector perp to both LOS and toroidal dir [channel][{x,y,z}] */ private double uPerp[][]; /** magnetic field on los in [channel][x,y,z][point] */ private double[][][] Bxyz; /** magnetic field on los in [channel][mag,theta,alpha][point] */ private double[][][] Bmta; /** electron density and temperature on los in [channel][point] */ private double[][] ne, Te; /** cold Omega with ne = 1 [channel][point][para,perp,perpperp] */ private double[][][] coldOmegaFact; /** addition to cold Omega to get warm omega, with ne=1 and tau=1 [channel][point][para,perp,perpperp] */ private double[][][] warmOmegaFact; final double uPhi[] = new double[]{ 0, 1, 0 }; private double neUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronDensity"); private double TeUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronTemperature"); public PolarimeterPredictionNode() { this("Polarimetry predictions"); } public PolarimeterPredictionNode(String name) { super(name); addConnectionPoint(new ConnectionPoint("magneticField", MagneticField.class, false, getField("magneticField"))); addConnectionPoint(new ConnectionPoint("electronDensity", ScalarND.class, false, getField("electronDensity"))); addConnectionPoint(new ConnectionPoint("electronTemperature", ScalarND.class, false, getField("electronTemperature"))); addConnectionPoint(new ConnectionPoint("los", LineOfSightSystem.class, false, getField("los"))); addConnectionPoint(new ConnectionPoint("dPsiStatus", IntegerArray.class, false, getField("dPsiStatus"))); addConnectionPoint(new ConnectionPoint("chiStatus", IntegerArray.class, false, getField("chiStatus"))); addConnectionPoint(new ConnectionPoint("psi0", DoubleArray.class, false, getField("psi0"))); addConnectionPoint(new ConnectionPoint("chi0", DoubleArray.class, false, getField("chi0"))); addConnectionPoint(new ConnectionPoint("laserWavelength", DoubleValue.class, false, getField("laserWavelength"))); addConnectionPoint(new ConnectionPoint("returnPass", BooleanArray.class, false, getField("returnPass"))); } /** Calcs LOS vectors uPara and uPerp */ private void updateLOS(){ losPoints = los.getPoints(); losDists = los.getDistances(); nChannels = los.getNumLOS(); doublePass = returnPass.getBooleanArray(); uPara = new double[nChannels][]; uPerp = new double[nChannels][]; for(int i=0; i < nChannels; i++){ /* one vector parallel to los */ uPara[i] = new double[] { losPoints[i][0][1] - losPoints[i][0][0], losPoints[i][1][1] - losPoints[i][1][0], losPoints[i][2][1] - losPoints[i][2][0] }; double dl = Math.sqrt((uPara[i][0] * uPara[i][0]) + (uPara[i][1] * uPara[i][1]) + (uPara[i][2] * uPara[i][2])); uPara[i][0] /= dl; uPara[i][1] /= dl; uPara[i][2] /= dl; /* and one perp to both los and toroidal dir: uPerp = (u x phi) / |u x phi| */ uPerp[i] = new double[] { -uPara[i][2], 0, uPara[i][0] }; dl = Math.sqrt((uPerp[i][0] * uPerp[i][0]) + (uPerp[i][2] * uPerp[i][2])); if(dl <= 0) throw new RuntimeException("LOS is parallel to toroidal dir (y) os polarisation is ill defined."); uPerp[i][0] /= dl; uPerp[i][2] /= dl; } } /** Resets the channel prediction enable flags to all automatic */ private void initChannelEnableRequests() { /* If it already exists and is the right size, leave as is */ if(predEnableRequest != null && predEnableRequest.length == nChannels) return; predEnableRequest = new int[nChannels]; for(int i=0; i < nChannels; i++) predEnableRequest[i] = Diagnostic.AUTOMATIC; } /** Recalculates calcDPsi/calcChi - whether or not to predict a channel based on the current * state of the predictionEnables and the data statuses dPsiStatus/chiStatus. */ private void recalcCalculationEnable() { int nChannels = los.getNumLOS(); calcDPsi = new boolean[nChannels]; calcChi = new boolean[nChannels]; int statusDPsi[] = dPsiStatus.getIntegerArray(); int statusChi[] = chiStatus.getIntegerArray(); for(int i=0; i < nChannels; i++){ calcDPsi[i] = ((predEnableRequest[i] == Diagnostic.AUTOMATIC) && (statusDPsi[i] == Diagnostic.DATA_VALID )) || ((predEnableRequest[i] == Diagnostic.ON)); calcChi[i] = ((predEnableRequest[i] == Diagnostic.AUTOMATIC) && (statusChi[i] == Diagnostic.DATA_VALID )) || ((predEnableRequest[i] == Diagnostic.ON)); } } /** Checks that psi0 and chi0 are valid for the analytic approximated calcuation * and sets psiIs45Degrees[] accordingly */ private void checkInitialPolarisation(){ double psi0[] = this.psi0.getDoubleArray(); double chi0[] = this.chi0.getDoubleArray(); psiIs45Degrees = new boolean[nChannels]; for(int i=0; i < nChannels; i++){ if( (chi0[i]*chi0[i]) > (5*5*d2r*d2r) ) throw new IllegalArgumentException("Analytic approx not valid for chi0 != 0"); /* Get psi0 in degrees in 0 < psi0 < 180 */ double psiDegrees = psi0[i]*r2d; while(psiDegrees < 0) psiDegrees += 180; while(psiDegrees > 180) psiDegrees -= 180; if( (psiDegrees > 40 && psiDegrees < 50) || (psiDegrees > 130 && psiDegrees < 140)) psiIs45Degrees[i] = true; else if( psiDegrees < 5 || psiDegrees > 175 || (psiDegrees > 85 && psiDegrees < 95)) /* psi0 ~= 0 */ psiIs45Degrees[i] = false; else{ System.err.println("WARNING: Analytic approx not valid for psi0 = "+psiDegrees+" != {0 or +/-45} on channel "+i+". Disabling this channel's prediction. "); calcDPsi[i] = false; calcChi[i] = false; } } } /** Predicts dPsi and chi from full integration. * ds.dz = W ^ s * * W = * */ private void doFullPrediction() { double psi0[] = this.psi0.getDoubleArray(); double chi0[] = this.chi0.getDoubleArray(); dPsi = new double[nChannels]; chi = new double[nChannels]; for(int i=0;i