0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 #include <math.h>
0042 #include "mex.h"
0043
0044
0045
0046 #define ec_IN prhs[0]
0047 #define e_IN prhs[1]
0048 #define t_IN prhs[2]
0049 #define Z_IN prhs[3]
0050
0051
0052
0053 #define seEH_OUT plhs[0]
0054 #define ec_OUT plhs[1]
0055 #define e_OUT plhs[2]
0056 #define t_OUT plhs[3]
0057
0058
0059 static xgamma(x,y,re,im)
0060 double x,y,*re,*im;
0061 {
0062
0063 double w,re2,im2,fn,xs,m,at;
0064 double xq,yq,zq,z4,z6,h,h1,ln,a;
0065 double rekl,imkl,r,s,res,ims,re1,im1;
0066
0067 w = 2.5066283;
0068
0069 if (fabs(x)-3.0 < 0) {
0070 if (x < 0) {
0071 xs = x-3.0;
0072 m = 1;
0073 } else if (x >= 0) {
0074 xs = x+3.0;
0075 m = 2;
0076 }
0077 } else if (fabs(x)-3.0 >= 0) {
0078 xs = x;
0079 m = 3;
0080 }
0081
0082 xq = xs*xs;
0083 yq = y*y;
0084 zq = xq+yq;
0085 z4 = zq*zq;
0086 z6 = zq*z4;
0087 h = xq-yq;
0088 h1 = 1.0+xs/(12.0*zq);
0089 ln = .5*log(zq);
0090 at = atan2(y,xs);
0091 rekl = h1+h/(288.0*z4)-2.6813272e-3*xs*(xq-3.0*yq)/z6-2.294721e-4*(h*h-4.0*xq*yq)/(z4*z4);
0092 imkl = -y/(12.0*zq)*(h1-3.2175926e-2*(h+xq+xq)/z4-1.10146605e-2*h*xs/z6);
0093 r = w*exp(ln*(xs-.5)-xs-y*at);
0094 s = y*ln-y+(xs-.5)*at;
0095 res = cos(s);
0096 ims = sin(s);
0097 re1 = r*(rekl*res-imkl*ims);
0098 im1 = r*(rekl*ims+imkl*res);
0099
0100 if (m < 2) {
0101 re2 = (xs+1.0)*(xq+2.0*xs-3.0*yq);
0102 im2 = y*(3.0*xq-yq+6.0*xs+2.0);
0103 re[0] = re1*re2-im1*im2;
0104 im[0] = re1*im2+re2*im1;
0105 } else if (m == 2) {
0106 re2 = xs*(xq-3.0*yq+11.0)-6.0*h-6.0;
0107 im2 = y*(3.0*xq-yq-12.0*xs+11.0);
0108 fn = pow(re2,2)+pow(im2,2);
0109 re[0] = (re1*re2+im1*im2)/fn;
0110 im[0] = (re2*im1-re1*im2)/fn;
0111 } else if (m > 2) {
0112 re[0] = re1;
0113 im[0] = im1;
0114 }
0115 }
0116
0117
0118
0119 static
0120 #ifdef __STDC__
0121 hgfn(
0122 double a,
0123 double b,
0124 double c,
0125 double d,
0126 double e,
0127 double f,
0128 double x,
0129 double *re,
0130 double *im
0131 )
0132 #else
0133 hgfn(a,b,c,d,e,f,x,re,im)
0134 double a,b,c,d,e,f,x,*re,*im;
0135 #endif
0136 {
0137 double delta,h1,h2,h3,h4,h5,h6,h7;
0138 double h8,h9,h10,h11,h12,h13;
0139
0140 delta = (double)1.0e-07;
0141
0142 if (fabs(x) < 1.0) {
0143 h1 = a;
0144 h2 = c;
0145 h3 = e;
0146 h4 = (double)1.0;
0147 h5 = -b*d;
0148 h6 = f*f;
0149 h7 = a*c+h5;
0150 h8 = a*d+b*c;
0151 h9 = x/(e*e+h6);
0152 h10 = (h7*e+h8*f)*h9;
0153 re[0] = 1.0+h10;
0154 h11 = (h8*e-h7*f)*h9;
0155 im[0] = h11;
0156 h1 = h1+1.0;
0157 h2 = h2+1.0;
0158 h3 = h3+1.0;
0159 h4 = h4+1.0;
0160 h7 = h1*h2+h5;
0161 h8 = h1*d+h2*b;
0162 h9 = x/((h3*h3+h6)*h4);
0163 h12 = (h7*h3+h8*f)*h9;
0164 h13 = (h8*h3-h7*f)*h9;
0165 h7 = h10*h12-h11*h13;
0166 re[0] = re[0]+h7;
0167 h11 = h11*h12+h10*h13;
0168 im[0] = im[0]+h11;
0169 h10 = h7;
0170
0171 while (fabs(h10) > delta || fabs(h11) > delta) {
0172 h1 = h1+1.0;
0173 h2 = h2+1.0;
0174 h3 = h3+1.0;
0175 h4 = h4+1.0;
0176 h7 = h1*h2+h5;
0177 h8 = h1*d+h2*b;
0178 h9 = x/((h3*h3+h6)*h4);
0179 h12 = (h7*h3+h8*f)*h9;
0180 h13 = (h8*h3-h7*f)*h9;
0181 h7 = h10*h12-h11*h13;
0182 re[0] = re[0]+h7;
0183 h11 = h11*h12+h10*h13;
0184 im[0] = im[0]+h11;
0185 h10 = h7;
0186 }
0187 }
0188 }
0189
0190
0191 static
0192 #ifdef __STDC__
0193 hgf(
0194 double a,
0195 double b,
0196 double c,
0197 double d,
0198 double g,
0199 double x,
0200 double *re,
0201 double *im
0202 )
0203 #else
0204 hgf(a,b,c,d,g,x,re,im)
0205 double a,b,c,d,g,x,*re,*im;
0206 #endif
0207 {
0208
0209 double delta,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11;
0210
0211 delta = (double)1.0e-07;
0212
0213 if(fabs(x) < 1.0) {
0214 h1 = a;
0215 h2 = c;
0216 h3 = g;
0217 h4 = (double)1.0;
0218 h5 = -b*d;
0219 h6 = (a*c+h5)*x/g;
0220 re[0] = (double)1.0+h6;
0221 h7 = (a*d+b*c)*x/g;
0222 im[0] = h7;
0223 h1 = h1+(double)1.0;
0224 h2 = h2+(double)1.0;
0225 h3 = h3+(double)1.0;
0226 h4 = h4+(double)1.0;
0227 h8 = x/(h3*h4);
0228 h9 = h1*h2+h5;
0229 h10 = h1*d+h2*b;
0230 h11 = (h6*h9-h7*h10)*h8;
0231 re[0] = re[0]+h11;
0232 h7 = (h7*h9+h6*h10)*h8;
0233 im[0] = im[0]+h7;
0234 h6 = h11;
0235
0236 while (fabs(h6) > delta || fabs(h7) > delta) {
0237 h1 = h1+(double)1.0;
0238 h2 = h2+(double)1.0;
0239 h3 = h3+(double)1.0;
0240 h4 = h4+(double)1.0;
0241 h8 = x/(h3*h4);
0242 h9 = h1*h2+h5;
0243 h10 = h1*d+h2*b;
0244 h11 = (h6*h9-h7*h10)*h8;
0245 re[0] = re[0]+h11;
0246 h7 = (h7*h9+h6*h10)*h8;
0247 im[0] = im[0]+h7;
0248 h6 = h11;
0249
0250 }
0251 }
0252 }
0253
0254
0255
0256 static
0257 #ifdef __STDC__
0258 sigma(
0259 double *z,
0260 double *t1,
0261 double *pen,
0262 double *theta1,
0263 double *sec
0264 )
0265 #else
0266 sigma(z,t1,pen,theta1,sec)
0267 double *z,*t1,*pen,*theta1,*sec;
0268 #endif
0269 {
0270 double alpha,mcq,r0q,g1,g2,pi,xbarn;
0271 double arc,gnte,dte,gnp,dph,e1,t2,e2,beta2,a;
0272 double eps1,eps1q,p1q,p1,eps2,eps2q,p2q;
0273 double p2,p12,k,kq,mu,rho,kappa,a1,a2,a21;
0274 double regn,gnq,reg,faktor;
0275 double th1,si1,cs1,p1k,d1,eta1q,w2q;
0276 double t,the,sie,cse,p1p2,p,phi,cs2,pp,p2k;
0277 double eta2q,eta12,d2,q2,ex,x;
0278 double al,r1,csr1,sir1,ex1,r2,csr2,sir2;
0279 double hi1,hi2,hi3,hi4,rev,imv,hi5,hi6,rew,imw;
0280 double rea1,ima1,rea2,ima2,reb,imb,ee1,ee2,e3,o1;
0281 double f1,f2,f3,wq,wq1,wq2,reg1,reg2,reg21;
0282 double ref1,ref2,ref3,re,im,wqint;
0283 double var0,var1,var2,var3,a21m,a1m,a2m,xex1m;
0284
0285 double zz,tt1,ppen,ttheta1;
0286
0287 int nte,np,it,ip,jt,jp,kk = 0;
0288
0289 double img1,img2,img21,imgn,img,imf1,imf2,imf3;
0290 double wqs,xn,wqskal;
0291
0292 alpha = (double)7.29736e-03;
0293 mcq = (double)511.006;
0294 r0q = (double)0.079408;
0295 g1 = (double)0.422649731;
0296 g2 = (double)1.577350269;
0297 pi = (double)3.14159265352;
0298 xbarn = (double)1.0e-24;
0299 arc = pi/(double)180.0;
0300 nte = 100;
0301 gnte = (double)100.0;
0302 dte = pi/gnte;
0303 np = 20;
0304 gnp = (double)20;
0305 dph = pi/gnp;
0306
0307
0308 e1 = t1[0];
0309 t2 = t1[0]-pen[0];
0310 e2 = t2;
0311
0312
0313 beta2 = (double)1.0-1.0/(pow((double)1.0+t1[0]/mcq,(double)2.0));
0314
0315 a = alpha*z[0];
0316
0317 eps1 = (double)1.0+e1/mcq;
0318 eps1q = pow(eps1,(double)2.0);
0319 p1q = eps1q-1.0;
0320 p1 = sqrt(p1q);
0321 eps2 = (double)1.0+e2/mcq;
0322 eps2q = pow(eps2,(double)2.0);
0323 p2q = eps2q-(double)1.0;
0324 p2 = sqrt(p2q);
0325 p12 = p1*p2;
0326 k = eps1-eps2;
0327 kq = k*k;
0328 mu = pow(p1+p2,(double)2.0)-kq;
0329 rho = (double)1.0/p1+(double)1.0/p2;
0330 kappa = eps1/p1+eps2/p2;
0331 a1 = a*eps1/p1;
0332 a2 = a*eps2/p2;
0333 a21 = a2-a1;
0334
0335 var0 = 0.0;
0336
0337 xgamma (var0,a1,®1,&img1);
0338 xgamma (var0,a2,®2,&img2);
0339 xgamma (var0,a21,®21,&img21);
0340
0341
0342 regn = reg1*reg2+img1*img2;
0343 imgn = reg1*img2-img1*reg2;
0344 gnq = (double)1.0/(pow(regn,(double)2.0)+pow(imgn,(double)2.0));
0345 reg = (reg21*regn+img21*imgn)*gnq;
0346 img = (regn*img21-imgn*reg21)*gnq;
0347 faktor = (double)4.0*a1*a2/((exp((double)2.0*pi*a1)-(double)1.0)*((double)1.0-exp(-2.0*pi*a2)))*r0q*a*z[0]*p2*k/p1;
0348
0349 th1 = theta1[0]*arc;
0350 si1 = sin(th1);
0351 cs1 = cos(th1);
0352 p1k = p1*k*cs1;
0353 d1 = (double)2.0*(eps1*k-p1k);
0354 eta1q = pow(p1*si1,(double)2.0);
0355 wq2 = (double)0.0;
0356
0357 for (it = 2; it<=nte; it += 2) {
0358 t = it*(double)1.0;
0359 the = dte*(t-g1);
0360 for (jt = 1; jt<=2; jt++) {
0361 sie = sin(the);
0362 cse = cos(the);
0363 p1p2 = p12*cse;
0364 wq1 = (double)0.0;
0365 for (ip = 2; ip<=np; ip += 2) {
0366 p=ip*(double)1.0;
0367 phi=dph*(p-g1);
0368 for (jp = 1; jp<=2; jp++) {
0369 cs2 = cs1*cse+si1*sie*cos(phi);
0370 pp = p12*cs1*cs2;
0371 p2k = p2*k*cs2;
0372 eta2q = p2q*(1.0-pow(cs2,(double)2.0));
0373 eta12 = p1p2-pp;
0374 d2 = (double)2.0*(eps2*k-p2k);
0375 q2 = d1+(double)2.0*(p2q-p1p2+p2k);
0376 ex = mu*q2/(d1*d2);
0377 x = (double)1.0-ex;
0378
0379 if ((x+(double)0.4) < 0.0) {
0380 al = log(ex);
0381 r1 = a1*al;
0382 csr1 = cos(r1);
0383 sir1 = sin(r1);
0384 ex1 = (double)1.0/ex;
0385 if((x+(double)1.5) < 0.0) {
0386 r2 = a2*al;
0387 csr2 = cos(r2);
0388 sir2 = sin(r2);
0389
0390 var0 = (double)0.0;
0391 var1 = (double)1.0;
0392 var2 = (double)1.0;
0393 a21m = -a21;
0394 a1m = -a1;
0395 a2m = -a2;
0396
0397 hgfn(var0,a1,var1,a2m,var2,a21m,ex1,&ref1,&imf1);
0398 hgfn(var0,a2,var1,a1m,var2,a21,ex1,&ref2,&imf2);
0399
0400
0401 hi1 = ref1*csr1+imf1*sir1;
0402 hi2 = ref2*csr2+imf2*sir2;
0403 hi3 = ref1*sir1-imf1*csr1;
0404 hi4 = ref2*sir2-imf2*csr2;
0405 rev = (reg*hi3-img*hi1)/a1+(reg*hi4+img*hi2)/a2;
0406 imv = (reg*hi1+img*hi3)/a1+(reg*hi2-img*hi4)/a2;
0407
0408 var0 = (double)1.0;
0409 var1 = (double)1.0;
0410 var2 = (double)1.0;
0411 a21m = -a21;
0412 a2m = -a2;
0413
0414 hgfn(var0,a1,var1,a2m,var2,a21m,ex1,&ref3,&imf3);
0415
0416 hi5 = reg*ref3-img*imf3;
0417 hi6 = reg*imf3+img*ref3;
0418 rew = (hi5*(csr1+csr2)+hi6*(sir1-sir2))/(a1*a2*ex);
0419 imw = (hi6*(csr1-csr2)-hi5*(sir1+sir2))/(a1*a2*ex);
0420
0421 rea1 = (rev+a1*ex*imw)/(d1*q2);
0422 ima1 = (imv-a1*ex*rew)/(d1*q2);
0423 rea2 = (rev+a2*ex*imw)/(d2*q2);
0424 ima2 = (imv-a2*ex*rew)/(d2*q2);
0425 reb = -a*imw/(d1*d2);
0426 imb = a*rew/(d1*d2);
0427 ee1 = ((double)4.0*eps2q-q2)*eta1q+((double)2.0*kq/d2*(eta2q+(double)1.0)+eta1q-eta12)*d1;
0428 ee2 = ((double)4.0*eps1q-q2)*eta2q+((double)2.0*kq/d1*(eta1q+(double)1.0)-eta2q+eta12)*d2;
0429 e3 = ((double)4.0*eps1*eps2-q2)*eta12+(double)0.5*d1*(eta12-eta2q)+0.5*d2*(eta1q-eta12)+2.0*kq*(eta12+(double)1.0);
0430 o1 = (double)0.5*mu*((double)1.0/p2-(double)1.0/p1)*pp;
0431 f1 = k*(rho*(eta1q-eta12)-o1+(double)0.25*mu*d1/p1+kappa*p1*cs1*p1p2+(2.0/p1-kappa*p12/k)*p2k+k*(eps2*rho+kappa));
0432 f2 = k*(rho*(eta2q-eta12)+o1-(double)0.25*mu*d2/p2+kappa*p2*cs2*p1p2-(2.0/p2+kappa*p12/k)*p1k-k*(eps1*rho+kappa));
0433 f3 = mu*(kq-kq*cs1*cs2+(p1q-p2q)/(p1q*p2q)*(p1q-p2q+p1k+p2k))-2.0*kq*rho*rho;
0434 wq = faktor*(ee1*(pow(rea1,(double)2.0)+pow(ima1,(double)2.0))+ee2*(pow(rea2,(double)2.0)+pow(ima2,(double)2.0))
0435 -(double)2.0*(e3*(rea1*rea2+ima1*ima2)+f1*(rea1*reb+ima1*imb)
0436 +f2*(rea2*reb+ima2*imb))+f3*(pow(reb,(double)2.0)+pow(imb,(double)2.0)));
0437 wq1 = wq1+wq;
0438
0439 }
0440
0441 else {
0442
0443 var0 = (double)0.0;
0444 var1 = (double)1.0;
0445 var2 = (double)1.0;
0446 a2m = -a2;
0447 xex1m = -x*ex1;
0448
0449 hgf(var0,a1,var1,a2m,var2,xex1m,&re,&im);
0450
0451 rev = re*csr1+im*sir1;
0452 imv = im*csr1-re*sir1;
0453
0454 var0 = (double)1.0;
0455 var1 = (double)1.0;
0456 var2 = (double)2.0;
0457 a2m = -a2;
0458 xex1m = -x*ex1;
0459
0460 hgf(var0,a1,var1,a2m,var2,xex1m,&re,&im);
0461
0462 rew = (re*csr1+im*sir1)*ex1;
0463 imw = (im*csr1-re*sir1)*ex1;
0464 rea1 = (rev+a1*ex*imw)/(d1*q2);
0465 ima1 = (imv-a1*ex*rew)/(d1*q2);
0466 rea2 = (rev+a2*ex*imw)/(d2*q2);
0467 ima2 = (imv-a2*ex*rew)/(d2*q2);
0468 reb = -a*imw/(d1*d2);
0469 imb = a*rew/(d1*d2);
0470 ee1 = ((double)4.0*eps2q-q2)*eta1q+(2.0*kq/d2*(eta2q+(double)1.0)+eta1q-eta12)*d1;
0471 ee2 = ((double)4.0*eps1q-q2)*eta2q+(2.0*kq/d1*(eta1q+(double)1.0)-eta2q+eta12)*d2;
0472 e3 = ((double)4.0*eps1*eps2-q2)*eta12+0.5*d1*(eta12-eta2q)+0.5*d2*(eta1q-eta12)+2.0*kq*(eta12+(double)1.0);
0473 o1 = (double)0.5*mu*((double)1.0/p2-(double)1.0/p1)*pp;
0474 f1 = k*(rho*(eta1q-eta12)-o1+(double)0.25*mu*d1/p1+kappa*p1*cs1*p1p2+((double)2.0/p1-kappa*p12/k)*p2k+k*(eps2*rho+kappa));
0475 f2 = k*(rho*(eta2q-eta12)+o1-(double)0.25*mu*d2/p2+kappa*p2*cs2*p1p2-((double)2.0/p2+kappa*p12/k)*p1k-k*(eps1*rho+kappa));
0476 f3 = mu*(kq-kq*cs1*cs2+(p1q-p2q)/(p1q*p2q)*(p1q-p2q+p1k+p2k))-(double)2.0*kq*rho*rho;
0477 wq = faktor*(ee1*(pow(rea1,(double)2.0)+pow(ima1,(double)2.0))+ee2*(pow(rea2,(double)2.0)+pow(ima2,(double)2.0))
0478 -(double)2.0*(e3*(rea1*rea2+ima1*ima2)+f1*(rea1*reb+ima1*imb)
0479 +f2*(rea2*reb+ima2*imb))+f3*(pow(reb,(double)2.0)+pow(imb,(double)2.0)));
0480 wq1 = wq1+wq;
0481 }
0482 }
0483 else {
0484
0485 var0 = (double)0.0;
0486 var1 = (double)0.0;
0487 var2 = (double)1.0;
0488
0489 hgf(var0,a1,var1,a2,var2,x,&rev,&imv);
0490
0491 var0 = (double)1.0;
0492 var1 = (double)1.0;
0493 var2 = (double)2.0;
0494
0495 hgf(var0,a1,var1,a2,var2,x,&rew,&imw);
0496
0497 rea1 = (rev+a1*ex*imw)/(d1*q2);
0498 ima1 = (imv-a1*ex*rew)/(d1*q2);
0499 rea2 = (rev+a2*ex*imw)/(d2*q2);
0500 ima2 = (imv-a2*ex*rew)/(d2*q2);
0501 reb = -a*imw/(d1*d2);
0502 imb = a*rew/(d1*d2);
0503 ee1 = ((double)4.0*eps2q-q2)*eta1q+(2.0*kq/d2*(eta2q+(double)1.0)+eta1q-eta12)*d1;
0504 ee2 = ((double)4.0*eps1q-q2)*eta2q+(2.0*kq/d1*(eta1q+(double)1.0)-eta2q+eta12)*d2;
0505 e3 = ((double)4.0*eps1*eps2-q2)*eta12+0.5*d1*(eta12-eta2q)+0.5*d2*(eta1q-eta12)+(double)2.0*kq*(eta12+(double)1.0);
0506 o1 = (double)0.5*mu*((double)1.0/p2-(double)1.0/p1)*pp;
0507 f1 = k*(rho*(eta1q-eta12)-o1+(double)0.25*mu*d1/p1+kappa*p1*cs1*p1p2+((double)2.0/p1-kappa*p12/k)*p2k+k*(eps2*rho+kappa));
0508 f2 = k*(rho*(eta2q-eta12)+o1-(double)0.25*mu*d2/p2+kappa*p2*cs2*p1p2-((double)2.0/p2+kappa*p12/k)*p1k-k*(eps1*rho+kappa));
0509 f3 = mu*(kq-kq*cs1*cs2+(p1q-p2q)/(p1q*p2q)*(p1q-p2q+p1k+p2k))-2.0*kq*rho*rho;
0510 wq = faktor*(ee1*(pow(rea1,(double)2.0)+pow(ima1,(double)2.0))+ee2*(pow(rea2,(double)2.0)+pow(ima2,(double)2.0))
0511 -(double)2.0*(e3*(rea1*rea2+ima1*ima2)+f1*(rea1*reb+ima1*imb)
0512 +f2*(rea2*reb+ima2*imb))+f3*(pow(reb,(double)2.0)+pow(imb,(double)2.0)));
0513 wq1 = wq1+wq;
0514 }
0515 phi=dph*(p-g2);
0516 }
0517 }
0518 wq2 = wq2+wq1*sie;
0519 the = dte*(t-g2);
0520 }
0521 }
0522 wqint = 2.0*wq2*dte*dph;
0523 wqs = wqint/mcq;
0524 wqskal = wqs*beta2*pen[0]/(z[0]*z[0]);
0525 xn = xbarn*wqs*sqrt(beta2);
0526 sec[0] = xn*mcq/sqrt(beta2);
0527 }
0528
0529
0530
0531
0532
0533 static
0534 #ifdef __STDC__
0535 sigmabrem(
0536 double seEHout[],
0537 double ecout[],
0538 double eout[],
0539 double tout[],
0540 double ec[],
0541 double e[],
0542 double t[],
0543 double *Z,
0544 unsigned int m,
0545 unsigned int n,
0546 unsigned int mode
0547 )
0548 #else
0549 sigmabrem(seEHout,ecout,eout,tout,ec,e,t,Z,m,n,mode)
0550 double seEHout[],ecout[],eout[],tout[],ec[],e[],t[],*Z;
0551 unsigned int m,n,mode;
0552 #endif
0553 {
0554
0555 unsigned int k,i,j;
0556 double pi = 3.14159265352;
0557 double sec,ee,ephot,ang;
0558
0559 if (mode == 1) {
0560
0561 k = 0;
0562 for (i = 0; ifor (j = 0; jif (ecout[k] <= eout[k])
0573 seEHout[k] = 0;
0574 else {
0575 sigma(Z,&ee,&ephot,&ang,&sec);
0576 seEHout[k] = sec;
0577 }
0578 }
0579 }
0580
0581 } else if (mode == 2) {
0582
0583 k = 0;
0584 for (i = 0; ifor (j = 0; jif (ecout[k] <= eout[k])
0594 seEHout[k] = 0;
0595 else {
0596 sigma(Z,&ee,&ephot,&ang,&sec);
0597 seEHout[k] = sec;
0598 }
0599
0600 }
0601 }
0602
0603 } else if (mode == 3) {
0604
0605 k = 0;
0606 for (i = 0; ifor (j = 0; jif (ecout[k] <= eout[k])
0616 seEHout[k] = 0;
0617 else {
0618 sigma(Z,&ee,&ephot,&ang,&sec);
0619 seEHout[k] = sec;
0620 }
0621 }
0622 }
0623
0624 }
0625 }
0626
0627
0628
0629
0630
0631 void mexFunction( int nlhs,
0632 mxArray *plhs[],
0633 int nrhs,
0634 const mxArray *prhs[]
0635 )
0636 {
0637
0638 double *seEHout,*ecout,*eout,*tout;
0639
0640 double *ec,*e,*t,*Z;
0641 unsigned int m,n,mec,nec,me,ne,mt,nt;
0642 unsigned int mode = 0;
0643
0644
0645
0646 if (nrhs < 2) {
0647
0648 mexErrMsgTxt("The function eh_dke_yp requires at least two input arguments.");
0649
0650 } else if (nrhs == 2) {
0651
0652 ec = mxGetPr(ec_IN);
0653 e = mxGetPr(e_IN);
0654
0655 t = (double *)mxCalloc(1,sizeof(double));
0656 Z = (double *)mxCalloc(1,sizeof(double));
0657
0658 *t = (double) 0.0;
0659 *Z = (double) 1.0;
0660
0661 mec = mxGetM(ec_IN);
0662 nec = mxGetN(ec_IN);
0663 me = mxGetM(e_IN);
0664 ne = mxGetN(e_IN);
0665 mt = 1;
0666 nt = 1;
0667
0668 } else if (nrhs == 3) {
0669
0670 ec = mxGetPr(ec_IN);
0671 e = mxGetPr(e_IN);
0672 t = mxGetPr(t_IN);
0673
0674 Z = (double *)mxCalloc(1,sizeof(double));
0675
0676 *Z = (double) 1.0;
0677
0678 mec = mxGetM(ec_IN);
0679 nec = mxGetN(ec_IN);
0680 me = mxGetM(e_IN);
0681 ne = mxGetN(e_IN);
0682 mt = mxGetM(t_IN);
0683 nt = mxGetN(t_IN);
0684
0685 } else if (nrhs == 4) {
0686
0687 ec = mxGetPr(ec_IN);
0688 e = mxGetPr(e_IN);
0689 t = mxGetPr(t_IN);
0690 Z = mxGetPr(Z_IN);
0691
0692 mec = mxGetM(ec_IN);
0693 nec = mxGetN(ec_IN);
0694 me = mxGetM(e_IN);
0695 ne = mxGetN(e_IN);
0696 mt = mxGetM(t_IN);
0697 nt = mxGetN(t_IN);
0698
0699 } else if (nlhs > 3) {
0700 mexErrMsgTxt("The function eh_dke_yp requires up to three output arguments.");
0701 }
0702
0703
0704
0705
0706
0707 if (fmin(mec,nec) != 1 || fmin(me,ne) != 1 || fmin(mt,nt) != 1) {
0708 mexErrMsgTxt("The function eh_dke_yp requires that ec and e be vectors.");
0709 }
0710
0711 if (fmax(mec,nec) != 1 && fmin(me,ne) != 1 && fmin(mt,nt) != 1) {
0712 mexErrMsgTxt("The function eh_dkeyp requires only two vectors as input arguments.");
0713 }
0714
0715 if (fmax(mt,nt) == 1) {
0716 mode = 1;
0717 n = fmax(mec,nec);
0718 m = fmax(me,ne);
0719 }
0720 if (fmax(me,ne) == 1) {
0721 mode = 2;
0722 n = fmax(mec,nec);
0723 m = fmax(mt,nt);
0724 }
0725 if (fmax(mec,nec) == 1) {
0726 mode = 3;
0727 n = fmax(mt,nt);
0728 m = fmax(me,ne);
0729 }
0730
0731
0732
0733
0734 seEH_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
0735 ec_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
0736 e_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
0737 t_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
0738
0739
0740
0741 seEHout = mxGetPr(seEH_OUT);
0742 ecout = mxGetPr(ec_OUT);
0743 eout = mxGetPr(e_OUT);
0744 tout = mxGetPr(t_OUT);
0745
0746
0747
0748 sigmabrem(seEHout,ecout,eout,tout,ec,e,t,Z,m,n,mode);
0749
0750
0751
0752 if (nrhs == 2) {
0753 mxFree(t);
0754 mxFree(Z);
0755 } else if (nrhs == 3)
0756 mxFree(Z);
0757
0758 return;
0759
0760 }
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782