Friday 25 February 2022

A Global Warming Calculation. Doubling CO2 to 800ppmv makes one fifth of a degree warming.

So since around COP26, in some of my spare time. I calculated just how much CO2 (and water vapour), absorbs heat radiation from earth. Based on the 3 (4 for H20) strongest absorption lines in the HITRAN database that come from absorption of a photon by the ground state. Going from 400 ppvm CO2 to 800 ppvm, increases absorbtion by just 1.1 Watts per square meter. Below is the chart of absorption by part per million watts of CO2 in the atmosphere.
I have made the code public, at my github. Finishing off a paper with the details of the calculation, but it is basically integrating over latitude, frequency and height, the absorption line and line width with a Lorentzian Line shape. Lorentz is known to be larger than actually absorption at frequencies away from the peak, so this is an over estimate. So how much heating does the extra 1.1 Watts cause. If all the heat goes into the emission, then by Stefan Boltzaman, https://en.wikipedia.org/wiki/Climate_sensitivity, Note 2, 1 Watt is 1 Kelvin (or C) warming. But that for any warming, some of that heat going into evaporating water, so the temperature change well be lower. I will look at evaporation at a later time. Updating this post, now I have done climate sensitivity including rainfall calculation, while the naive temperature increase is about 1.1Kelvin, with rainfall it is only a sixth of that about 0.167 Kelvin.

Wednesday 23 February 2022

100 days of data at Dresden-II on elastic neutrino interaction with nucleii.

https://arxiv.org/abs/2202.10829, Bounds on new physics with data of the Dresden-II reactor experiment and COHERENT Coloma et al.
Therefore, current data from these experiments is not enough by itself to impose meaningful constraints on the complete parameter space of NSI with quarks, if considered in full generality
Maybe in 3 years, a thousand days of data, will yield something on non standard interaction. The above paper does yield limits on neutrino magnetic moments. Another paper today,https://arxiv.org/abs/2202.10622, Implications of the first evidence for coherent elastic scattering of reactor neutrinos, Liao, Lui and Marfatia. Yields.
If the standard Lindhard model correctly describes the quenching factor, the data may indicate a light vector or scalar mediator, or a large neutrino magnetic moment. We await more data

Thursday 10 February 2022

Need a Multithreaded Numerical Integrator in Java - See below

In Java, the quickest Thread Exector for Thread that run many times in the Fork Join Pool. You can speed up your numerical integration on multicored CPU like most computers are today. Here it is in Git Hub If we define a function
public abstract class DoubFunction {


    abstract double evalInner(double x, double params[], int i);

    public double eval(double x, int i, double ...params){
        return evalInner(x,  params, i);
    }

}

Then example Simpson Rule Integrates in Java are. With 16 Java Threads on AMD Ryzen with 8 Cores 16 CPU Threads, Standard Integrate takes 51.4 sec, and the ThreadIntegrator only 32.5 seconds.
import java.util.Arrays;
import java.util.List;
import java.util.concurrent.*;
import java.util.concurrent.atomic.AtomicReference;

public class SimpsonsRule {

    private static double third = 1.0/3.0;

    private static int THREADS = 16;

    public static double integrateThreaded(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        ForkJoinPool pool = new ForkJoinPool(THREADS);
       AtomicDouble sum = new AtomicDouble();
        for(int i=0; i {
                double mul = ii%2==0? 2*third: 4*third;
                if (ii==0) { mul = third; }
                if (ii==N-1){ mul = third; }
                double mul1 = mul;
                double x = a + h * ii;
                double fi = f.eval(x,ii, params);
                if (Double.isNaN(fi) ){
                    System.err.println(f.getClass().getName() + "IS NaN at "+x);
                }
                sum.getAndAdd(mul1*fi);
            });
        }
        try {
            pool.shutdown();
            if (!pool.awaitTermination(1, TimeUnit.HOURS)){
                pool.shutdownNow();
            }
        } catch (Exception e){}
        return sum.getAndAdd(0.0) * h;
    }

    public static double integrate(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        double fa = f.eval(a,0, params);
        double fb = f.eval(b, N-1, params);
        if (Double.isNaN(fa) ){
          System.err.println(f.getClass().getName() + "IS NaN at "+a);
        }
        if (Double.isNaN(fb)){
            System.err.println(f.getClass().getName() + "IS NaN at "+b);
        }
        // 1/3 terms
        double sum = third * (fa + fb);


        // 4/3 terms
        for (int i = 1; i < N - 1; i += 2) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            sum += 4.0 * third * fx;
        }

        // 2/3 terms
        for (int i = 2; i < N - 1; i += 2) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            sum += 2.0 * third * fx;
        }

        return sum * h;
    }

    public static double integrateConsecutive(double a, double b, int N, DoubFunction f, double ...params) {         // precision parameter
        double h = (b - a) / (N - 1);     // step size
        double fa = f.eval(a,0, params);
        if (Double.isNaN(fa) ){
            System.err.println(f.getClass().getName() + "IS NaN at "+a);
        }

        // 1/3 terms
        double sum = third*fa;


        // 4/3 terms
        boolean isOdd = true;
        for (int i = 1; i < N - 1; i += 1) {
            double x = a + h * i;
            double fx = f.eval(x,i, params);
            if (Double.isNaN(fx)){
                System.err.println(f.getClass().getName() + "IS NaN at "+x);
            }
            if (isOdd) {
                sum += 4.0 * third * fx;
            } else {
                sum += 2.0 * third * fx;
            }
            isOdd=!isOdd;
        }


        double fb = f.eval(b, N-1, params);
        if (Double.isNaN(fb)){
            System.err.println(f.getClass().getName() + "IS NaN at "+b);
        }
        sum = sum+ fb*third;

        return sum * h;
    }

    public static void main(String args[]){
        // Roots of polynumerial to integrate
        List in = Arrays.asList(-0.9, -0.8,-0.7, -0.6,-0.5, -0.4, -0.3, -0.2, -0.1, 0, .1,.2,.3,.4,.5, .6,.7, .8, .9 );
        DoubFunction func = new DoubFunction() {
            @Override
            double evalInner(double x, double[] params, int i) {
                return in.stream().map(y->y.doubleValue()-x).reduce(1.0,(a,b)->(a*b));
            }
        };
        double consec=0;
        long startConsec = System.currentTimeMillis();
        for(int i=1;i<1000; i++) {
            consec = integrateConsecutive(-1, 1, 100000, func);
        }
        double timeConsec = (System.currentTimeMillis() - startConsec)/1000.0;
        double standard=0;
        long startStandard = System.currentTimeMillis();
        for(int i=1; i<1000; i++) {
                standard = integrate(-1,1,100000,func);
        }
        double timeStandard = (System.currentTimeMillis() - startStandard)/1000.0;
        double threaded=0;
        long startThreaded = System.currentTimeMillis();
        for(int i=1;i<1000; i++) {
            threaded = integrateThreaded(-1, 1, 100000, func);
        }
        double timeThreaded = (System.currentTimeMillis() - startThreaded)/1000.0;
        System.out.println("Standard Integrator: "+standard+" time taken: "+timeStandard+" seconds");
        System.out.println("Consecutive Integrator: "+consec+" time taken: "+timeConsec+" seconds");
        System.out.println("threaded Integrator: "+threaded+" time taken: "+timeThreaded+" seconds");
    }

}

class AtomicDouble {
    private AtomicReference value = new AtomicReference(Double.valueOf(0.0));
    double getAndAdd(double delta) {
        while (true) {
            Double currentValue = value.get();
            Double newValue = Double.valueOf(currentValue.doubleValue() + delta);
            if (value.compareAndSet(currentValue, newValue))
                return currentValue.doubleValue();
        }
    }
}