//2 component vector to hold the real and imaginary parts of a complex number: typedef float2 cfloat; #define I ((cfloat)(0.0, 1.0)) /* * Return Real (Imaginary) component of complex number: */ inline float real(cfloat a){ return a.x; } inline float imag(cfloat a){ return a.y; } /* * Get the modulus of a complex number (its length): */ inline double cmod(cfloat a){ return (sqrt(a.x*a.x + a.y*a.y)); } /* * Get the argument of a complex number (its angle): */ /* inline float carg(cfloat a) { return atan(a.y/a.x); }*/ inline float carg(cfloat a){ if(a.x > 0){ return atan(a.y / a.x); }else if(a.x < 0 && a.y >= 0){ return atan(a.y / a.x) + M_PI_F; }else if(a.x < 0 && a.y < 0){ return atan(a.y / a.x) - M_PI_F; }else if(a.x == 0 && a.y > 0){ return -M_PI_F/2; }else if(a.x == 0 && a.y < 0){ return M_PI_F/2; }else{ return 0; } } /* * Multiply two complex numbers: * * a = (aReal + I*aImag) * b = (bReal + I*bImag) * a * b = (aReal + I*aImag) * (bReal + I*bImag) * = aReal*bReal +I*aReal*bImag +I*aImag*bReal +I^2*aImag*bImag * = (aReal*bReal - aImag*bImag) + I*(aReal*bImag + aImag*bReal) */ inline cfloat cmult(cfloat a, cfloat b){ return (cfloat)( a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x); } /*Square of one complex number: * * a = (aReal + I*aImag) * a^2=(aReal + I*aImag)*(aReal + I*aImag) * =(aReal*aReal-aImag*aImag)+I*(2*aReal*aImag) */ inline cfloat csqr(cfloat a){ return (cfloat)(a.x*a.x-a.y*a.y, 2*a.x*a.y); } /* * Divide two complex numbers: * * aReal + I*aImag (aReal + I*aImag) * (bReal - I*bImag) * ----------------- = --------------------------------------- * bReal + I*bImag (bReal + I*bImag) * (bReal - I*bImag) * * aReal*bReal - I*aReal*bImag + I*aImag*bReal - I^2*aImag*bImag * = --------------------------------------------------------------- * bReal^2 - I*bReal*bImag + I*bImag*bReal -I^2*bImag^2 * * aReal*bReal + aImag*bImag aImag*bReal - Real*bImag * = ---------------------------- + I* -------------------------- * bReal^2 + bImag^2 bReal^2 + bImag^2 * */ inline cfloat cdiv(cfloat a, cfloat b){ return (cfloat)((a.x*b.x + a.y*b.y)/(b.x*b.x + b.y*b.y), (a.y*b.x - a.x*b.y)/(b.x*b.x + b.y*b.y)); } /* * Square root of complex number. */ inline cfloat csqrt(cfloat a){ return (cfloat)( sqrt(cmod(a)) * cos(carg(a)/2), sqrt(cmod(a)) * fabs(sin(carg(a)/2))); } /*inline cfloat csqrt(cfloat a){ return (cfloat)( sqrt((cmod(a)+a.x)/2), sign(a.y)*sqrt((cmod(a)-a.x)/2)); }*/ /*Cosine and Sinus of complex number * *cos( aReal + I*aImag )=cos(aReal)*cosh(aImag)-I*sin(aReal)*sinh(aImag) *sin( aReal + I*aImag )=sin(aReal)*cosh(aImag)+I*cos(aReal)*sinh(aImag) * */ inline cfloat ccos(cfloat a){ return (cfloat) (cos(a.x)*cosh(a.y), -sin(a.x)*sinh(a.y)); } inline cfloat csin(cfloat a){ return (cfloat) (sin(a.x)*cosh(a.y), cos(a.x)*sinh(a.y)); } inline cfloat csinh(cfloat a){ return (cfloat) (cos(-a.y)*(exp(a.x)-exp(-a.x))/2,-sin(-a.y)*(exp(a.x)+exp(-a.x))/2); } /*Exponential function of complex number * *exp(aReal + I*aImag)=exp(aReal)*(cos(aImag)+I*sin(aImag)) * */ inline cfloat cexp(cfloat a){ return (cfloat) (exp(a.x)*cos(a.y),exp(a.x)*sin(a.y)); } kernel void XRR9381132d834748d8ba3821d6f1e060fd( global read_only double* params, global read_only double* X, global double* Y) { int index = get_global_id(0); double zcoordinate_0 = params[0*4+1]; double delta_0 = params[0*4+2]; double beta_0 = params[0*4+3]; double roughness_0 = params[0*4+4]; double zcoordinate_1 = params[1*4+1]; double delta_1 = params[1*4+2]; double beta_1 = params[1*4+3]; double roughness_1 = params[1*4+4]; double zcoordinate_2 = params[2*4+1]; double delta_2 = params[2*4+2]; double beta_2 = params[2*4+3]; double roughness_2 = params[2*4+4]; double zcoordinate_3 = params[3*4+1]; double delta_3 = params[3*4+2]; double beta_3 = params[3*4+3]; double roughness_3 = params[3*4+4]; const double K_squared=(2*M_PI_F/1.5418)*(2*M_PI_F/1.5418); const cfloat q_squared=csqr((cfloat)(X[index]/2,0)); cfloat K_0=csqrt(q_squared+(cfloat)(-2*K_squared*delta_0,2*K_squared*beta_0)); cfloat K_1=csqrt(q_squared+(cfloat)(-2*K_squared*delta_1,2*K_squared*beta_1)); cfloat K_2=csqrt(q_squared+(cfloat)(-2*K_squared*delta_2,2*K_squared*beta_2)); cfloat K_3=csqrt(q_squared+(cfloat)(-2*K_squared*delta_3,2*K_squared*beta_3)); cfloat FresnelCoeff_0 = cdiv((K_0 - K_1) , (K_0 + K_1)); cfloat FresnelCoeff_1 = cdiv((K_1 - K_2) , (K_1 + K_2)); cfloat FresnelCoeff_2 = cdiv((K_2 - K_3) , (K_2 + K_3)); cfloat tmp=csqr((cfloat)(M_PI_F,0))-(cfloat)(8,0); cfloat FunctionDivision_1 =cexp(cmult((cfloat)(-4*roughness_1*roughness_1,0),cmult(K_0,K_1))); cfloat FunctionDivision_2 =cexp(cmult((cfloat)(-4*roughness_2*roughness_2,0),cmult(K_1,K_2))); cfloat FunctionDivision_3 =cexp(cmult((cfloat)(-4*roughness_3*roughness_3,0),cmult(K_2,K_3))); cfloat ModifiedFresnelC_0 = cmult(FunctionDivision_1 , -FresnelCoeff_0); cfloat ModifiedFresnelC_1 = cmult(FunctionDivision_2 , -FresnelCoeff_1); cfloat ModifiedFresnelC_2 = cmult(FunctionDivision_3 , -FresnelCoeff_2); cfloat Ratio_3=(0,0); cfloat exp_minus_2=cexp(cmult((cfloat)(0,-2*(zcoordinate_2)),K_2)); cfloat exp_plus_2=cexp(cmult((cfloat)(0,2*(zcoordinate_2)),K_3)); cfloat Ratio_2=cmult(exp_minus_2,cdiv(ModifiedFresnelC_2+cmult(Ratio_3,exp_plus_2),(cfloat)(1,0)+cmult(cmult(ModifiedFresnelC_2,Ratio_3),exp_plus_2))); cfloat exp_minus_1=cexp(cmult((cfloat)(0,-2*(zcoordinate_1)),K_1)); cfloat exp_plus_1=cexp(cmult((cfloat)(0,2*(zcoordinate_1)),K_2)); cfloat Ratio_1=cmult(exp_minus_1,cdiv(ModifiedFresnelC_1+cmult(Ratio_2,exp_plus_1),(cfloat)(1,0)+cmult(cmult(ModifiedFresnelC_1,Ratio_2),exp_plus_1))); cfloat exp_minus_0=cexp(cmult((cfloat)(0,-2*(zcoordinate_0)),K_0)); cfloat exp_plus_0=cexp(cmult((cfloat)(0,2*(zcoordinate_0)),K_1)); cfloat Ratio_0=cmult(exp_minus_0,cdiv(ModifiedFresnelC_0+cmult(Ratio_1,exp_plus_0),(cfloat)(1,0)+cmult(cmult(ModifiedFresnelC_0,Ratio_1),exp_plus_0))); Y[index] =Ratio_0.x*Ratio_0.x+Ratio_0.y*Ratio_0.y; }