ok, here it comes again

. thanks for your answer!
_________________________________________________________
#include "mex.h"
#include <math.h>
#define PI 3.14159
void source (double x, double y, double z, double x1, double y1, double z1,
double x2, double y2, double z2, double x3, double y3, double
z3,
double x4, double y4, double z4, double *srcphi, double
*dblphi)
{
double eps = 1e-10;
double d12,d23,d34,d41;
double r1,r2,r3,r4;
double e1,e2,e3,e4;
double h1,h2,h3,h4;
double m12,m23,m34,m41;
x1=x1+1*eps;
x2=x2-1*eps;
x3=x3+2*eps;
x4=x4-2*eps;
z1=0;
z2=0;
z3=0;
z4=0;
d12 = sqrt(pow((x2-x1),2) + pow((y2-y1),2));
d23 = sqrt(pow((x3-x2),2) + pow((y3-y2),2));
d34 = sqrt(pow((x4-x3),2) + pow((y4-y3),2));
d41 = sqrt(pow((x1-x4),2) + pow((y1-y4),2));
r1 = sqrt(pow((x-x1),2) + pow((y-y1),2) + pow(z,2));
r2 = sqrt(pow((x-x2),2) + pow((y-y2),2) + pow(z,2));
r3 = sqrt(pow((x-x3),2) + pow((y-y3),2) + pow(z,2));
r4 = sqrt(pow((x-x4),2) + pow((y-y4),2) + pow(z,2));
e1 = pow((x-x1),2) + pow(z,2);
e2 = pow((x-x2),2) + pow(z,2);
e3 = pow((x-x3),2) + pow(z,2);
e4 = pow((x-x4),2) + pow(z,2);
h1 = (x-x1)*(y-y1);
h2 = (x-x2)*(y-y2);
h3 = (x-x3)*(y-y3);
h4 = (x-x4)*(y-y4);
m12 = (y2-y1)/(x2-x1);
m23 = (y3-y2)/(x3-x2);
m34 = (y4-y3)/(x4-x3);
m41 = (y1-y4)/(x1-x4);
*(srcphi) = -1.0/(4*PI)*(
(
(((x-x1)*(y2-y1)-(y-y1)*(x2-x1))/d12)*log((r1+r2+d12)/(r1+r2-d12)) +
(((x-x2)*(y3-y2)-(y-y2)*(x3-x2))/d23)*log((r2+r3+d23)/(r2+r3-d23)) +
(((x-x3)*(y4-y3)-(y-y3)*(x4-x3))/d34)*log((r3+r4+d34)/(r3+r4-d34)) +
(((x-x4)*(y1-y4)-(y-y4)*(x1-x4))/d41)*log((r4+r1+d41)/(r4+r1-d41))
)
+fabs(z)*(
atan((m12*e1-h1)/(z*r1)) - atan((m12*e2-h2)/(z*r2)) +
atan((m23*e2-h2)/(z*r2)) - atan((m23*e3-h3)/(z*r3)) +
atan((m34*e3-h3)/(z*r3)) - atan((m34*e4-h4)/(z*r4)) +
atan((m41*e4-h4)/(z*r4)) - atan((m41*e1-h1)/(z*r1))
)
);
*(dblphi) = -1.0/(4*PI)*(
atan((m12*e1-h1)/(z*r1)) - atan((m12*e2-h2)/(z*r2))
+ atan((m23*e2-h2)/(z*r2)) - atan((m23*e3-h3)/(z*r3))
+ atan((m34*e3-h3)/(z*r3)) - atan((m34*e4-h4)/(z*r4))
+ atan((m41*e4-h4)/(z*r4)) - atan((m41*e1-h1)/(z*r1)));
return;
}
void coefficientmatrix(int totpanelsf, int totpanelswake, double eps,
double *xcf, double *ycf, double *zcf,
double *x1f, double *y1f, double *z1f,
double *x2f, double *y2f, double *z2f,
double *x3f, double *y3f, double *z3f,
double *x4f, double *y4f, double *z4f,
double *xe1, double *ye1, double *ze1,
double *xe2, double *ye2, double *ze2,
double *nxf, double *nyf, double *nzf,
double *z1fqs, double *z2fqs, double *z3fqs, double
*z4fqs,
double *xcwake, double *ycwake, double *zcwake,
double *x1wake, double *y1wake, double *z1wake,
double *x2wake, double *y2wake, double *z2wake,
double *x3wake, double *y3wake, double *z3wake,
double *x4wake, double *y4wake, double *z4wake,
double *xe1wake, double *ye1wake, double *ze1wake,
double *xe2wake, double *ye2wake, double *ze2wake,
double *xe3wake, double *ye3wake, double *ze3wake,
double *z1wakeqs, double *z2wakeqs, double *z3wakeqs,
double *z4wakeqs,
double *A, double *S, double *E)
{
int i,j,k,count = 0,count2=0;
double xp, yp, zp;
double pcolloc[3], p1[3], p2[3], p3[3], p4[3];
double srcphi, dblphi;
for(i=0;i<totpanelsf;i++)
{
xp = xcf
-nxf*eps;
yp = ycf-nyf*eps;
zp = zcf-nzf*eps;
for (j=0;j<totpanelsf;j++)
{
pcolloc[0] = xe1[j]*(xp-xcf[j])+ye1[j]*(yp-ycf[j])+ze1[j]*(zp-zcf[j]);
pcolloc[1] = xe2[j]*(xp-xcf[j])+ye2[j]*(yp-ycf[j])+ze2[j]*(zp-zcf[j]);
pcolloc[2] = nxf[j]*(xp-xcf[j])+nyf[j]*(yp-ycf[j])+nzf[j]*(zp-zcf[j]);
p1[0] =
xe1[j]*(x1f[j]-xcf[j])+ye1[j]*(y1f[j]-ycf[j])+ze1[j]*(z1fqs[j]-zcf[j]);
p1[1] =
xe2[j]*(x1f[j]-xcf[j])+ye2[j]*(y1f[j]-ycf[j])+ze2[j]*(z1fqs[j]-zcf[j]);
p1[2] =
nxf[j]*(x1f[j]-xcf[j])+nyf[j]*(y1f[j]-ycf[j])+nzf[j]*(z1fqs[j]-zcf[j]);
p2[0] =
xe1[j]*(x2f[j]-xcf[j])+ye1[j]*(y2f[j]-ycf[j])+ze1[j]*(z2fqs[j]-zcf[j]);
p2[1] =
xe2[j]*(x2f[j]-xcf[j])+ye2[j]*(y2f[j]-ycf[j])+ze2[j]*(z2fqs[j]-zcf[j]);
p2[2] =
nxf[j]*(x2f[j]-xcf[j])+nyf[j]*(y2f[j]-ycf[j])+nzf[j]*(z2fqs[j]-zcf[j]);
p3[0] =
xe1[j]*(x3f[j]-xcf[j])+ye1[j]*(y3f[j]-ycf[j])+ze1[j]*(z3fqs[j]-zcf[j]);
p3[1] =
xe2[j]*(x3f[j]-xcf[j])+ye2[j]*(y3f[j]-ycf[j])+ze2[j]*(z3fqs[j]-zcf[j]);
p3[2] =
nxf[j]*(x3f[j]-xcf[j])+nyf[j]*(y3f[j]-ycf[j])+nzf[j]*(z3fqs[j]-zcf[j]);
p4[0] =
xe1[j]*(x4f[j]-xcf[j])+ye1[j]*(y4f[j]-ycf[j])+ze1[j]*(z4fqs[j]-zcf[j]);
p4[1] =
xe2[j]*(x4f[j]-xcf[j])+ye2[j]*(y4f[j]-ycf[j])+ze2[j]*(z4fqs[j]-zcf[j]);
p4[2] =
nxf[j]*(x4f[j]-xcf[j])+nyf[j]*(y4f[j]-ycf[j])+nzf[j]*(z4fqs[j]-zcf[j]);
source(pcolloc[0],pcolloc[1],pcolloc[2],p1[0],p1[1],p1[2],
p2[0],p2[1],p2[2],p3[0],p3[1],p3[2],
p4[0],p4[1],p4[2],&srcphi, &dblphi);
*(S+count)= srcphi;
*(A+count)= dblphi;
count++;
}
for (k=0;k<totpanelswake;k++)
{
pcolloc[0] =
xe1wake[k]*(xp-xcwake[k])+ye1wake[k]*(yp-ycwake[k])+ze1wake[k]*(zp-zcwake[k]
);
pcolloc[1] =
xe2wake[k]*(xp-xcwake[k])+ye2wake[k]*(yp-ycwake[k])+ze2wake[k]*(zp-zcwake[k]
);
pcolloc[2] =
xe3wake[k]*(xp-xcwake[k])+ye3wake[k]*(yp-ycwake[k])+ze3wake[k]*(zp-zcwake[k]
);
p1[0] =
xe1wake[k]*(x1wake[k]-xcwake[k])+ye1wake[k]*(y1wake[k]-ycwake[k])+ze1wake[k]
*(z1wakeqs[k]-zcwake[k]);
p1[1] =
xe2wake[k]*(x1wake[k]-xcwake[k])+ye2wake[k]*(y1wake[k]-ycwake[k])+ze2wake[k]
*(z1wakeqs[k]-zcwake[k]);
p1[2] =
xe3wake[k]*(x1wake[k]-xcwake[k])+ye3wake[k]*(y1wake[k]-ycwake[k])+ze3wake[k]
*(z1wakeqs[k]-zcwake[k]);
p2[0] =
xe1wake[k]*(x2wake[k]-xcwake[k])+ye1wake[k]*(y2wake[k]-ycwake[k])+ze1wake[k]
*(z2wakeqs[k]-zcwake[k]);
p2[1] =
xe2wake[k]*(x2wake[k]-xcwake[k])+ye2wake[k]*(y2wake[k]-ycwake[k])+ze2wake[k]
*(z2wakeqs[k]-zcwake[k]);
p2[2] =
xe3wake[k]*(x2wake[k]-xcwake[k])+ye3wake[k]*(y2wake[k]-ycwake[k])+ze3wake[k]
*(z2wakeqs[k]-zcwake[k]);
p3[0] =
xe1wake[k]*(x3wake[k]-xcwake[k])+ye1wake[k]*(y3wake[k]-ycwake[k])+ze1wake[k]
*(z3wakeqs[k]-zcwake[k]);
p3[1] =
xe2wake[k]*(x3wake[k]-xcwake[k])+ye2wake[k]*(y3wake[k]-ycwake[k])+ze2wake[k]
*(z3wakeqs[k]-zcwake[k]);
p3[2] =
xe3wake[k]*(x3wake[k]-xcwake[k])+ye3wake[k]*(y3wake[k]-ycwake[k])+ze3wake[k]
*(z3wakeqs[k]-zcwake[k]);
p4[0] =
xe1wake[k]*(x4wake[k]-xcwake[k])+ye1wake[k]*(y4wake[k]-ycwake[k])+ze1wake[k]
*(z4wakeqs[k]-zcwake[k]);
p4[1] =
xe2wake[k]*(x4wake[k]-xcwake[k])+ye2wake[k]*(y4wake[k]-ycwake[k])+ze2wake[k]
*(z4wakeqs[k]-zcwake[k]);
p4[2] =
xe3wake[k]*(x4wake[k]-xcwake[k])+ye3wake[k]*(y4wake[k]-ycwake[k])+ze3wake[k]
*(z4wakeqs[k]-zcwake[k]);
source(pcolloc[0],pcolloc[1],pcolloc[2],p1[0],p1[1],p1[2],
p2[0],p2[1],p2[2],p3[0],p3[1],p3[2],
p4[0],p4[1],p4[2],&srcphi, &dblphi);
*(E+count2)= dblphi;
count2++;
}
}
return;
}
void mexFunction (int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
/* Input variables */
int totpanelsf, totpanelswake;
double eps;
double *xcf, *ycf, *zcf;
double *x1f, *y1f, *z1f;
double *x2f, *y2f, *z2f;
double *x3f, *y3f, *z3f;
double *x4f, *y4f, *z4f;
double *xe1, *ye1, *ze1;
double *xe2, *ye2, *ze2;
double *nxf, *nyf, *nzf;
double *z1fqs, *z2fqs, *z3fqs, *z4fqs;
double *xcwake, *ycwake, *zcwake;
double *x1wake, *y1wake, *z1wake;
double *x2wake, *y2wake, *z2wake;
double *x3wake, *y3wake, *z3wake;
double *x4wake, *y4wake, *z4wake;
double *xe1wake, *ye1wake, *ze1wake;
double *xe2wake, *ye2wake, *ze2wake;
double *xe3wake, *ye3wake, *ze3wake;
double *z1wakeqs, *z2wakeqs, *z3wakeqs, *z4wakeqs;
/* Output variables */
double *A, *S, *E;
/* Assign scalars to input */
totpanelsf = mxGetScalar(prhs[0]);
totpanelswake = mxGetScalar(prhs[1]);
eps = mxGetScalar(prhs[2]);
/* Assign pointers to input */
xcf = mxGetPr(prhs[3]);
ycf = mxGetPr(prhs[4]);
zcf = mxGetPr(prhs[5]);
x1f = mxGetPr(prhs[6]);
y1f = mxGetPr(prhs[7]);
z1f = mxGetPr(prhs[8]);
x2f = mxGetPr(prhs[9]);
y2f = mxGetPr(prhs[10]);
z2f = mxGetPr(prhs[11]);
x3f = mxGetPr(prhs[12]);
y3f = mxGetPr(prhs[13]);
z3f = mxGetPr(prhs[14]);
x4f = mxGetPr(prhs[15]);
y4f = mxGetPr(prhs[16]);
z4f = mxGetPr(prhs[17]);
xe1 = mxGetPr(prhs[18]);
ye1 = mxGetPr(prhs[19]);
ze1 = mxGetPr(prhs[20]);
xe2 = mxGetPr(prhs[21]);
ye2 = mxGetPr(prhs[22]);
ze2 = mxGetPr(prhs[23]);
nxf = mxGetPr(prhs[24]);
nyf = mxGetPr(prhs[25]);
nzf = mxGetPr(prhs[26]);
z1fqs = mxGetPr(prhs[27]);
z2fqs = mxGetPr(prhs[28]);
z3fqs = mxGetPr(prhs[29]);
z4fqs = mxGetPr(prhs[30]);
xcwake = mxGetPr(prhs[31]);
ycwake = mxGetPr(prhs[32]);
zcwake = mxGetPr(prhs[33]);
x1wake = mxGetPr(prhs[34]);
y1wake = mxGetPr(prhs[35]);
z1wake = mxGetPr(prhs[36]);
x2wake = mxGetPr(prhs[37]);
y2wake = mxGetPr(prhs[38]);
z2wake = mxGetPr(prhs[39]);
x3wake = mxGetPr(prhs[40]);
y3wake = mxGetPr(prhs[41]);
z3wake = mxGetPr(prhs[42]);
x4wake = mxGetPr(prhs[43]);
y4wake = mxGetPr(prhs[44]);
z4wake = mxGetPr(prhs[45]);
xe1wake = mxGetPr(prhs[46]);
ye1wake = mxGetPr(prhs[47]);
ze1wake = mxGetPr(prhs[48]);
xe2wake = mxGetPr(prhs[49]);
ye2wake = mxGetPr(prhs[50]);
ze2wake = mxGetPr(prhs[51]);
xe3wake = mxGetPr(prhs[52]);
ye3wake = mxGetPr(prhs[53]);
ze3wake = mxGetPr(prhs[54]);
z1wakeqs = mxGetPr(prhs[55]);
z2wakeqs = mxGetPr(prhs[56]);
z3wakeqs = mxGetPr(prhs[57]);
z4wakeqs = mxGetPr(prhs[58]);
/* Create Matrices for the outputs */
plhs[0] = mxCreateDoubleMatrix(totpanelsf,totpanelsf,mxREAL);
plhs[1] = mxCreateDoubleMatrix(totpanelsf,totpanelsf,mxREAL);
plhs[2] = mxCreateDoubleMatrix(totpanelsf,totpanelswake,mxREAL);
/* Assign pointers to output */
A = mxGetPr(plhs[0]);
S = mxGetPr(plhs[1]);
E = mxGetPr(plhs[2]);
/* Call the coefficient matrix C subroutine */
coefficientmatrix(totpanelsf,totpanelswake,eps,
xcf, ycf, zcf,
x1f, y1f, z1f,
x2f, y2f, z2f,
x3f, y3f, z3f,
x4f, y4f, z4f,
xe1, ye1, ze1,
xe2, ye2, ze2,
nxf, nyf, nzf,
z1fqs, z2fqs, z3fqs, z4fqs,
xcwake, ycwake, zcwake,
x1wake, y1wake, z1wake,
x2wake, y2wake, z2wake,
x3wake, y3wake, z3wake,
x4wake, y4wake, z4wake,
xe1wake, ye1wake, ze1wake,
xe2wake, ye2wake, ze2wake,
xe3wake, ye3wake, ze3wake,
z1wakeqs, z2wakeqs, z3wakeqs, z4wakeqs,
A, S, E);
return;
}