volume = abs(1.0f/6.0f*((X[1]-X[0])*((Y[1]-Y[2])*(Z[2]-Z[3])-(Y[2]-Y[3])*(Z[1]-Z[2])) + (X[2]-X[1])*((Y[2]-Y[3])*(Z[0]-Z[1])-(Y[0]-Y[1])*(Z[2]-Z[3])) + (X[3]-X[2])*((Y[0]-Y[1])*(Z[1]-Z[2])-(Y[1]-Y[2])*(Z[0]-Z[1]))));

        for(int j = 0; j<4; ++j){
            ii = j;
            jj = (ii + 1)*(ii <= 2);
            kk = (jj + 1)*(jj<=2);
            ll = (kk + 1)*(kk<=2);

            Ai<< X[jj], Y[jj], Z[jj],
                X[kk], Y[kk], Z[kk],
                X[ll], Y[ll], Z[ll];

            Bi<< 1, Y[jj], Z[jj],
                1, Y[kk], Z[kk],
                1, Y[ll], Z[ll];

            Ci<< Y[jj], 1,  Z[jj],
                Y[kk], 1,  Z[kk],
                Y[ll], 1,  Z[ll];

            Di<< Y[jj], Z[jj], 1,
                Y[kk], Z[kk], 1,
                Y[ll], Z[ll], 1;

            a[ii] = Ai.determinant();
            b[ii] = -1*Bi.determinant();
            c[ii] = Ci.determinant();
            d[ii] = -1*Di.determinant();

            std::cout<<"b "<<b[ii]<<" ratio "<<1.0f/(2.0f*volume)<<" prod "<<(float)b[ii]*1.0f/(2.0f*volume)<<"\n";