import std.stdio, std.math,std.complex, std.conv,std.parallelism,std.algorithm,std.range;

alias Creal = Complex!real;

///fysical constants in our problem
immutable static real sigma=56_000_000.0;
immutable static real vacPermOver2Pi=2*0.000_000_1; //permeability = 4*pi*10^-7
immutable static real copperRelPerm= 0.999994;

///the configuration of our problem. Changing this implies recalculating L
immutable static real width=0.002;
immutable static real height=0.002;
immutable static real distance=0.001;

///a magical number that analytically cancels out but can potentially boost precision by eliminating small numbers
immutable static real compensator=1.0;

///the function we have to integrate
real analyticalresult(real x1,real x2,real a){
    a=(a==0)?0.0000000000001:a;//cheap trick to avoid numerical problems
    return 1.0/4*(3*x1*x1 - 6*x1*x2 + x2*x2 - 4*a*(x1-x2)*atan((x1-x2)/a) + (a*a - (x1-x2)*(x1-x2))*log(a*a+(x1-x2)*(x1-x2)));
}

void main()
{
    writeln("generating L things");
    Grid problemGrid=simpleGrid(3,3,width,height,distance);
    auto N = problemGrid.length;
    auto L = generateLmatrix(problemGrid);


    foreach(freq;[100,1_000_000,10_000_000,100_000_000,1_000_000_000,10_000_000_000,100_000_000_000]){
        writeln("generating data for freq = ",freq);

        auto omega = 2*PI*freq.to!real();

        auto coefMatrix = new Creal[][](N+1,N+1);
        coefMatrix[N][N]=Creal(0.0,0.0);

        for(int row=0;row<N;row++)
            for(int col=0;col<N;col++)
                coefMatrix[row][col]=Creal(0.0,-omega*vacPermOver2Pi*L[row][col]);

        for(int i=0;i<N;i++){
            coefMatrix[i][i]+=Creal(problemGrid[i].area()*compensator/sigma); //diagonal
            coefMatrix[N][i]=Creal(1.0*problemGrid[i].area()*compensator); //lowest row
            coefMatrix[i][N]=Creal(-100.0*problemGrid[i].area()*compensator); //furthest right col
        }

        auto inv=computeInverse(coefMatrix);

        Creal[] currents;
        for(int i=0;i<N;i++)
            currents~=inv[i][N]*100.0*compensator;//*problemGrid[i].area();
        //saveCurrentVisualization(problemGrid,currents,freq.to!string()); //not in slimem

        auto concoction = inv[N][N]*100.0*compensator; //this is R+jwL
        writeln("R : ",concoction.re * 1_000," miliohm/meter");
        writeln("L : ",concoction.im/omega*1_000_000_000, " nanohenri/meter");

        writeln();

    }

    writeln("press enter you twat");
    stdin.readln();
}

real[][] generateLmatrix(const ref Grid problemGrid){
    auto N = problemGrid.length;

    real[][] L = new real[][](N,N);

    foreach(row/*index*/,ref real[] currow;taskPool.parallel(L)){
        PointTuple p1=problemGrid[row];
        foreach(ind,ref curcel;currow){
            PointTuple p2=problemGrid[ind];
            //readability is only useful if you plan  to read it
            curcel=integrate2d( (y1,y2) => compensator*(analyticalresult(p1.endx,p2.endx,abs(y1-y2)) - analyticalresult(p1.endx,p2.startx,abs(y1-y2)) -analyticalresult(p1.startx,p2.endx,abs(y1-y2)) + analyticalresult(p1.startx,p2.startx,abs(y1-y2))
                                                            -1*(analyticalresult(p1.endx,p2.endx,abs(y1+y2)) - analyticalresult(p1.endx,p2.startx,abs(y1+y2)) -analyticalresult(p1.startx,p2.endx,abs(y1+y2)) + analyticalresult(p1.startx,p2.startx,abs(y1+y2))))
                        , p1.starty, p1.endy, p2.starty, p2.endy);
            writeln("(",row,",",ind,") : ",curcel);
        }
    }
    return L;
}

///calculates the L21 norm of the given matrix abs(T element) has to be defined
real getNorm(T)(T[][] mat){
    real retval=0.0;
    for(int col=0;col<mat.length;col++){
        real rowThing=0.0;
        for(int row=0;row<mat.length;row++){
            rowThing+=abs(mat[row][col])*abs(mat[row][col]);
        }
        retval+=sqrt(rowThing);
    }
    return retval;
}

///using partial pivoting
///Todo : templatize? (not that trivial)
Creal [][] computeInverse(Creal [][] mat){
    auto inverse=new Creal [][] (mat.length,mat.length);
    for(int i=0;i<mat.length;i++)
        for(int j=0;j<mat.length;j++)
            inverse[i][j]=(i==j)?Creal(1.0):Creal(0.0);

    for(int i=0;i<mat.length;i++){
        //search the biggest spiel
        auto newSpiel=Creal(0.0);
        auto newSpielIndex=0;
        for(int iterator=i;iterator<mat.length;iterator++){
            if(abs(mat[iterator][i])>abs(newSpiel)){
                newSpiel=mat[iterator][i];
                newSpielIndex=iterator;
            }
        }

        if(abs(newSpiel)==0.0) // should be 'almost zero'
            throw new Exception("you dun goofed up");

        //rotate our biggest spiel-row to the right place
        for(int col=0;col<mat.length;col++){
            auto matDup=mat[newSpielIndex][col];
            auto invDup=inverse[newSpielIndex][col];

            mat[newSpielIndex][col]=mat[i][col];
            inverse[newSpielIndex][col]=inverse[i][col];

            mat[i][col]=matDup;
            inverse[i][col]=invDup;
        }

        auto multfact=Creal(1.0)/newSpiel;
        for(int col=0;col<mat.length;col++){//make diag element 1
            mat[i][col]*=multfact;
            inverse[i][col]*=multfact;
        }
        for(int row=0;row<mat.length;row++){//make elements at col i zero above our spiel
            if(row==i)
                continue;
            multfact=mat[row][i];
            for(int col=0;col<mat.length;col++){
                mat[row][col]-=multfact*mat[i][col];
                inverse[row][col]-=multfact*inverse[i][col];
            }
        }

    }
    return inverse;
}

///verry naive implementation of an integration algorithm (riemann sums)
real integrate2d(in real delegate(real,real) fun,in real x1,in real x2,in real y1,in real y2,in int numpieces=100 /*actually the numpieces is this squared*/){
    assert(x1<x2);
    assert(y1<y2);

    real xSegmentLength=(x2-x1)/numpieces;
    real ySegmentLength=(y2-y1)/numpieces;

    real retval=0.0;
    for(int i=0;i<numpieces*numpieces;i++)
        retval+=fun(x1 + (i%numpieces)*xSegmentLength + xSegmentLength/2,y1 + (i/numpieces)*ySegmentLength + ySegmentLength/2);

    return retval* xSegmentLength*ySegmentLength; //by only multiplying by the area(<<1) at the end we can potentially save a little precision
}

struct PointTuple{
    real startx;
    real starty;
    real endx;
    real endy;
    real area (){return (endy-starty)*(endx-startx);}
}

///by abstracting everything away, we can create really exotic grids
alias PointTuple[] Grid;

Grid simpleGrid(in uint widthPieces,in uint heightPieces,in real width,in real height,in real distance){
    Grid retval;
    for(int number=0;number<widthPieces*heightPieces;number++){
        PointTuple toret= PointTuple(width/widthPieces*(number/heightPieces),height/heightPieces*(number%heightPieces)+distance);
        toret.endx=toret.startx+width/widthPieces;
        toret.endy=toret.starty+height/heightPieces;
        retval~=toret;
    }
    return retval;
}