d2cross

double d2cross(double ***m, int ii, int jj, int kk, int dir1,int dir2,  int sh, int sm)
/*             matrix     , point                 , directions of dervs,shift , sample */
/*                                                , 1, 2, 3           , 0-left, 3,5,7 */
{         //order ==0 (first),dx=dx[dir1]*dx[dir2]
double tmp=0.0;
int i1,i2;
int dirr=6-dir1-dir2;
switch (sm*dirr) {
    case 3 : for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii][jj+i1-sh][kk+i2-sh]*kf3[0][sh][i1]*kf3[0][sh][i2]; break;
    case 6 : for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj][kk+i2-sh]*kf3[0][sh][i1]*kf3[0][sh][i2]; break;
    case 9 : for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj+i2-sh][kk]*kf3[0][sh][i1]*kf3[0][sh][i2]; break;
    case 5 : for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii][jj+i1-sh][kk+i2-sh]*kf5[0][sh][i1]*kf5[0][sh][i2]; break;
    case 10: for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj][kk+i2-sh]*kf5[0][sh][i1]*kf5[0][sh][i2]; break;
    case 15: for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj+i2-sh][kk]*kf5[0][sh][i1]*kf5[0][sh][i2]; break;
    case 7 : for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii][jj+i1-sh][kk+i2-sh]*kf7[0][sh][i1]*kf7[0][sh][i2]; break;
    case 14: for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj][kk+i2-sh]*kf7[0][sh][i1]*kf7[0][sh][i2]; break;
    case 21: for(i1=0; i1<sm; i1++)
             for(i2=0; i2<sm; i2++)
               tmp += m[ii+i1-sh][jj+i2-sh][kk]*kf7[0][sh][i1]*kf7[0][sh][i2]; break;
    default :
    nrerror("\nNO SUCH SAMPLE for derivative. Bye ...",0,0);
    }
return(tmp/dx[dir1-1]/dx[dir2-1]);
}

использование:

dA11[l][0][1] = dA11[l][1][0] = d2cross(f[l],i,j,k,2,1,(approx-1)/2,approx);

описание переменной:

double dA11[7][3][3];

описание для head.h:
double d2cross(double ***m,int ii,int jj,int kk,int dir1,int dir2,int sh,int sm);