/* * To change this template, choose Tools | Templates * and open the template in the editor. */ package org.das2.qds.demos; import java.text.ParseException; import java.util.Random; import java.util.logging.Level; import java.util.logging.Logger; import org.das2.datum.Datum; import org.das2.datum.DatumUtil; import org.das2.datum.DatumVector; import org.das2.datum.Units; import org.das2.qds.math.PoissonDistribution; import org.das2.qds.DDataSet; import org.das2.qds.QDataSet; import org.das2.qds.util.DataSetBuilder; /** * Model of plasma distribution function for given density, temperature, speed. * A java.util.Random object is passed in so that the data may be reproducible * (by using a given starting seed). * @author jbf */ public class PlasmaModel { private double nc = 1.2; // core density 1/cm^3 private double wcperp = 8000. * 1e5; // core thermal velocity cm/s private double wcparl = 8000. * 1e5; //static final double k = 1.380649e-16; // erg / k //static final double mass = 511000; // MeV //9.11e-28; // gram //static final double k = 8.617333262145e-5; // eV / k private static final double mass = 9.11e-28; // gram private double geomFactor = 1000e-40; PoissonDistribution poissonDistribution; public PlasmaModel() { poissonDistribution= new PoissonDistribution(); } /** * set the model density * @param density */ public void setDensity(Datum density) { this.nc = density.doubleValue(Units.pcm3); } /** * get the model density * @return */ public Datum getDensity() { return Units.pcm3.createDatum(nc); } /** * set the parallel speed * @param wcparl */ public void setWcparl(Datum wcparl) { this.wcparl = wcparl.doubleValue(Units.cmps); } /** * get the parallel speed. * @return */ public Datum getWcParl() { return Units.cmps.createDatum(this.wcparl); } /** * set the perpendicular speed * @param wcperp perpendicular speed */ public void setWcPerp(Datum wcperp) { this.wcperp= wcperp.doubleValue(Units.cmps); } /** * get the perpendicular speed * @return the perpendicular speed */ public Datum getWcPerp() { return Units.cmps.createDatum(this.wcperp); } /** * set the detector geometry factor * @param geom the detector geometry factor */ public void setGeomFactor( Datum geom ) { this.geomFactor= geom.doubleValue(Units.dimensionless); } /** * get the detector geometry factor * @return the detector geometry factor */ public Datum getGeomFactor( ) { return Units.dimensionless.createDatum(this.geomFactor); } /** * return f at the given energy * @param energy * @return */ public double f( Datum energy ) { double en= energy.doubleValue( Units.eV ); if ( wcperp!=wcparl ) { throw new IllegalArgumentException("distribution is not isotropic, need pitch angle"); } double v = Math.sqrt(2 * en * 1.6e-19 * 1e7 / mass); double logfc = Math.log10(nc / (Math.pow(Math.PI, (3. / 2.)) * wcparl * wcparl)) - 3 * Math.pow(v / wcparl, 2); return Math.pow(10,logfc); } public double f( Datum energy, Datum pitchAngle ) { double v = Math.sqrt(2 * energy.doubleValue(Units.eV) * 1.6e-19 * 1e7 / mass); double a = pitchAngle.doubleValue( Units.radians ); double vparl= Math.cos(a) * v; double vperp= Math.sin(a) * v; // double mass_2K= mass / ( 2 * k ); // double f= nc * ( mass_2K / ( Math.PI * wcperp ) ) // * Math.sqrt( mass_2K / ( Math.PI * wcperp ) ) // * Math.exp( - ( mass_2K * Math.pow( vperp, 2 ) / wcperp ) ) // * Math.exp( - ( mass_2K * Math.pow( vparl, 2 ) / wcparl ) ); //return f; double logfc = Math.log10(nc / (Math.pow(Math.PI, (3. / 2.)) * wcparl * wcperp)) - Math.pow(vparl / wcparl, 2) - 2*Math.pow(vperp / wcperp, 2); return Math.pow(10,logfc); } /** * return the counts at this energy, assuming an isotropic distribution. No * Poisson noise is added to the output. * @param energy * @return */ public double fcounts( Datum energy) { double en= energy.doubleValue(Units.eV); double fcount = 2. * (en / mass) * (en / mass) * geomFactor * f(energy); return fcount; } /** * return the counts at this energy, assuming an isotropic distribution, and * Poisson noise is added to the result. * @param energy in eV * @param random source of random numbers. * @return */ public int counts( Datum energy, Random random) { double en= energy.doubleValue(Units.eV); double fcount = 2. * (en / mass) * (en / mass) * geomFactor * f(energy); return poissonDistribution.poisson(fcount, random); } /** * return the counts at this energy and pitch angle, without Poisson noise added to the result. * @param energy * @param pitch * @return the floating point count rate. */ public double fcounts( Datum energy, Datum pitch ) { double f= f(energy, pitch); double fcount = f * 2. * Math.pow( energy.doubleValue( Units.eV) / mass, 2 ) * geomFactor; return fcount; } /** * return the counts at this energy and pitch angle, with Poisson noise added to the result. * @param energy * @param pitch * @param random random number source * @return the count rate */ public int counts(Datum energy, Datum pitch, Random random) { double fcount = fcounts( energy, pitch ); return poissonDistribution.poisson(fcount, random); } /** * return a rank 2 dataset with time as DEPEND_0 and energy as DEPEND_1. * @return */ public QDataSet getRank2( ) { try { Random random = new Random(5330); //System.err.println("First Random: "+random.nextDouble() ); //TODO: look into bug where this always causes hang in autoplot-test148 //System.err.println("Java version: "+System.getProperty("java.version")); Units xunits = Units.us2000; Datum start = Units.us2000.parse("2000-017T00:00"); Datum end = Units.us2000.parse("2000-018T00:00"); double xTagWidth = Units.seconds.convertDoubleTo( Units.microseconds, 13.8 ); double x = start.doubleValue(xunits); DataSetBuilder builder = new DataSetBuilder(2,1000,20); DataSetBuilder xx= new DataSetBuilder(1,1000); boolean ylog = false; DatumVector[] yTags = new DatumVector[1]; Random s = new Random(234567); // repeatable random sequence double n= 2.0; //int irec=0; while (x < end.doubleValue(xunits)) { int whichYTags = s.nextInt(yTags.length); int nj; if (yTags[whichYTags] == null) { nj = whichYTags * 10 + 20; double[] yy = new double[nj]; for (int j = 0; j < nj; j++) { if (ylog) { yy[j] = (nj / 300) + j * 0.05; // findbugs okay ICAST_IDIV_CAST_TO_DOUBLE yy[j] = Math.pow(10,yy[j]); } else { yy[j] = (nj / 3) + j * 1.2; } } yTags[whichYTags] = DatumVector.newDatumVector(yy, Units.eV); } else { nj = yTags[whichYTags].getLength(); } double[] zz = new double[nj]; int ncol = s.nextInt(4) + 1; DatumVector ydv = yTags[whichYTags]; for (int icol = 0; icol < ncol; icol++) { double d= random.nextDouble(); //if ( icol==0 && irec<150 ) { // System.err.println( String.format( "%d %d %5.4f", irun, irec, d ) ); //TODO: see above use at line 80. //} n= n * Math.pow(10, ( d-0.5 )/100 ); setDensity( Units.pcm3.createDatum(n) ); for (int j = 0; j < nj; j++) { zz[j] = counts( ydv.get(j), random ); builder.putValue( -1, j, zz[j] ); } xx.putValue( -1, x ); x += xTagWidth; //irec++; } builder.nextRecord(); xx.nextRecord(); } xx.putProperty( QDataSet.UNITS, Units.us2000 ); builder.putProperty( QDataSet.DEPEND_0, xx.getDataSet() ); DDataSet yy= DDataSet.wrap(yTags[0].toDoubleArray(Units.eV)); yy.putProperty(QDataSet.UNITS, Units.eV); builder.putProperty( QDataSet.DEPEND_1, yy ); //System.err.println("Last Random: "+random.nextDouble() ); //TODO: see above use at line 80. return builder.getDataSet(); } catch (ParseException ex) { throw new RuntimeException(ex); } } /** * return a rank 2 dataset with time as DEPEND_0 and energy as DEPEND_1. * @return */ public QDataSet getRank3( ) { try { Random random = new Random(5330); //System.err.println("First Random: "+random.nextDouble() ); //TODO: look into bug where this always causes hang in autoplot-test148 //System.err.println("Java version: "+System.getProperty("java.version")); Datum start = Units.us2000.parse("2000-017T00:00"); Datum end = Units.us2000.parse("2000-018T00:00"); Datum xTagWidth = Units.seconds.createDatum(13.8); Datum t = start; DataSetBuilder builder = new DataSetBuilder(3,1000,20,18); DataSetBuilder xx= new DataSetBuilder(1,1000); boolean ylog = false; DatumVector[] energyTags = new DatumVector[1]; DatumVector[] alphaTags = new DatumVector[1]; Random s = new Random(234567); // repeatable random sequence double n= 2.0; //int irec=0; while ( t.lt(end) ) { int whichTags = s.nextInt(energyTags.length); int ne; if (energyTags[whichTags] == null) { ne = whichTags * 10 + 20; double[] en = new double[ne]; for (int j = 0; j < ne; j++) { if (ylog) { en[j] = (ne / 300) + j * 0.05; // findbugs okay ICAST_IDIV_CAST_TO_DOUBLE en[j] = Math.pow(10,en[j]); } else { en[j] = (ne / 3) + j * 1.2; } } energyTags[whichTags] = DatumVector.newDatumVector(en, Units.eV); } else { ne = energyTags[whichTags].getLength(); } int na; if (alphaTags[whichTags] == null) { na = 18; double[] yy = new double[na]; for (int j = 0; j < na; j++) { yy[j]= Math.PI * ( 0.5 + j ) / na ; } alphaTags[whichTags] = DatumVector.newDatumVector(yy, Units.radians ); } else { na = alphaTags[whichTags].getLength(); } int ncol = s.nextInt(4) + 1; DatumVector ydv = energyTags[whichTags]; DatumVector adv = alphaTags[whichTags]; for (int icol = 0; icol < ncol; icol++) { double d= random.nextDouble(); //if ( icol==0 && irec<150 ) { // System.err.println( String.format( "%d %d %5.4f", irun, irec, d ) ); //TODO: see above use at line 80. //} n= n * Math.pow(10, ( d-0.5 )/100 ); setDensity( Units.pcm3.createDatum(n) ); for (int j = 0; j < ne; j++) { for ( int k=0; k