Révision 201 TrouNoir/TrouNoir.py
| TrouNoir.py (revision 201) | ||
|---|---|---|
| 86 | 86 |
float h,float r,float m,float bss) |
| 87 | 87 |
{
|
| 88 | 88 |
float flx; |
| 89 |
int fi; |
|
| 90 | 89 |
|
| 91 |
fi=(int)(rf*nbr/bss); |
|
| 92 |
flx=pow(r/m,q)*pow(rf,4)*b*db*h; |
|
| 90 |
flx=pow(r/m,q)*pow(rf,4)*b; |
|
| 93 | 91 |
return(flx); |
| 94 | 92 |
} |
| 95 | 93 |
|
| 96 | 94 |
float spectre_cn(float rf,float b,float db, |
| 97 | 95 |
float h,float r,float m,float bss) |
| 98 | 96 |
{
|
| 99 |
float flx;
|
|
| 100 |
float nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k;
|
|
| 97 |
double flx;
|
|
| 98 |
double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k,psc2,psk;
|
|
| 101 | 99 |
int fi,posfreq; |
| 102 | 100 |
|
| 103 | 101 |
planck=6.62e-34; |
| 104 | 102 |
k=1.38e-23; |
| 103 |
psk=5.38e-11; |
|
| 104 |
psc2=7.35e-51; |
|
| 105 | 105 |
temp=3.e7; |
| 106 | 106 |
m_point=1.e14; |
| 107 | 107 |
v=1.-3./r; |
| 108 | 108 |
|
| 109 |
qu=1/sqrt(1-3./r)/sqrt(r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))*(sqrt(6.)-sqrt(3.))/(sqrt(6.)+sqrt(3.))));
|
|
| 109 |
qu=1./sqrt(1.-3./r)/sqrt(r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))*(sqrt(6.)-sqrt(3.))/(sqrt(6.)+sqrt(3.))));
|
|
| 110 | 110 |
|
| 111 |
temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))* |
|
| 112 |
exp(-0.125*log(v))*exp(0.25*log(qu)); |
|
| 111 |
temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*exp(-0.125*log(v))*exp(0.25*log(qu)); |
|
| 113 | 112 |
|
| 114 | 113 |
flux_int=0; |
| 115 | 114 |
flx=0; |
| 116 | 115 |
|
| 117 |
for (fi=1;fi<nbr;fi++)
|
|
| 116 |
for (fi=0;fi<nbr;fi++)
|
|
| 118 | 117 |
{
|
| 119 |
nu_em=bss*fi/nbr;
|
|
| 118 |
nu_em=bss*(float)fi/(float)nbr;
|
|
| 120 | 119 |
nu_rec=nu_em*rf; |
| 121 |
posfreq=1./bss*nu_rec*nbr;
|
|
| 120 |
posfreq=(int)(1./bss*nu_rec*(float)nbr);
|
|
| 122 | 121 |
if ((posfreq>0)&&(posfreq<nbr)) |
| 123 | 122 |
{
|
| 124 |
flux_int=2*planck/9.e16f*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h; |
|
| 123 |
//flux_int=2*planck/9.e16f*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h; |
|
| 124 |
//flux_int=pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b; |
|
| 125 |
flux_int=2.*psc2*pow(nu_em,3)/(exp(psk*nu_em/temp_em)-1.)*pow(rf,3)*b; |
|
| 125 | 126 |
flx+=flux_int; |
| 126 | 127 |
} |
| 127 | 128 |
} |
| 129 |
//return(pow(rf,3)*b*temp_em); |
|
| 130 |
return((float)flx); |
|
| 128 | 131 |
return(flx); |
| 129 | 132 |
} |
| 130 | 133 |
|
| ... | ... | |
| 151 | 154 |
|
| 152 | 155 |
} |
| 153 | 156 |
|
| 154 |
__kernel void EachPixel(__global float *MyImage)
|
|
| 157 |
__kernel void EachPixel(__global float *zImage,__global float *fImage)
|
|
| 155 | 158 |
{
|
| 156 | 159 |
uint xi=(uint)get_global_id(0); |
| 157 | 160 |
uint yi=(uint)get_global_id(1); |
| ... | ... | |
| 185 | 188 |
//tho=PI/180.*80; |
| 186 | 189 |
tho=PI/180.*80.; |
| 187 | 190 |
q=-2; |
| 188 |
dist=1;
|
|
| 191 |
dist=0;
|
|
| 189 | 192 |
raie=0; |
| 190 | 193 |
|
| 191 | 194 |
nmx=(float)sizex/2.; |
| ... | ... | |
| 196 | 199 |
thx=asin(bmx/ro); |
| 197 | 200 |
pc=0; |
| 198 | 201 |
|
| 202 |
if (raie==0) |
|
| 203 |
{
|
|
| 204 |
bss=2.; |
|
| 205 |
} |
|
| 206 |
else |
|
| 207 |
{
|
|
| 208 |
bss=1.e6; |
|
| 209 |
} |
|
| 210 |
|
|
| 199 | 211 |
h=PI/5.e2f; |
| 200 | 212 |
|
| 201 | 213 |
// set origin as center of image |
| 202 | 214 |
float x=(float)xi-(float)(sizex/2)+(float)5e-1f; |
| 203 |
float y=(float)(sizey/2)-(float)yi-(float)5e-1f;
|
|
| 215 |
float y=(float)yi-(float)(sizey/2)+(float)5e-1f;
|
|
| 204 | 216 |
// angle extracted from cylindric symmetry |
| 205 | 217 |
phi=atanp(x,y); |
| 206 | 218 |
phd=atanp(cos(phi)*sin(tho),cos(tho)); |
| ... | ... | |
| 239 | 251 |
vp=vs; |
| 240 | 252 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
| 241 | 253 |
rp[nh]=fabs(b/us); |
| 242 |
// ExitOnImpact = (((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))||((fmod(ps,PI)<fmod(phd,PI))&&(fmod(pp,PI)>fmod(phd,PI))))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0; |
|
| 243 | 254 |
ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0; |
| 244 | 255 |
|
| 245 | 256 |
} while ((rp[nh]>=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0)); |
| ... | ... | |
| 250 | 261 |
else |
| 251 | 262 |
{
|
| 252 | 263 |
zp=0.; |
| 264 |
fp=0.; |
|
| 253 | 265 |
} |
| 254 | 266 |
|
| 255 | 267 |
// |
| 256 |
MyImage[yi+sizex*xi]=(float)zp; |
|
| 268 |
zImage[yi+sizex*xi]=(float)zp; |
|
| 269 |
fImage[yi+sizex*xi]=(float)fp; |
|
| 257 | 270 |
|
| 258 |
barrier(CLK_LOCAL_MEM_FENCE);
|
|
| 271 |
barrier(CLK_GLOBAL_MEM_FENCE);
|
|
| 259 | 272 |
|
| 260 | 273 |
} |
| 261 | 274 |
""") |
| ... | ... | |
| 269 | 282 |
image = Image.fromarray(SigmaInt) |
| 270 | 283 |
image.save("%s.jpg" % prefix)
|
| 271 | 284 |
|
| 272 |
def BlackHole(MyImage,Device):
|
|
| 285 |
def BlackHole(zImage,fImage,Device):
|
|
| 273 | 286 |
|
| 274 | 287 |
kernel_params = {}
|
| 275 | 288 |
|
| ... | ... | |
| 297 | 310 |
# Je recupere les flag possibles pour les buffers |
| 298 | 311 |
mf = cl.mem_flags |
| 299 | 312 |
|
| 300 |
print(MyImage)
|
|
| 313 |
print(zImage)
|
|
| 301 | 314 |
|
| 302 | 315 |
# Program based on Kernel2 |
| 303 |
MyImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=MyImage) |
|
| 316 |
zImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=zImage) |
|
| 317 |
fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) |
|
| 304 | 318 |
|
| 305 | 319 |
start_time=time.time() |
| 306 | 320 |
|
| 307 |
CLLaunch=BlackHoleCL.EachPixel(queue,(MyImage.shape[0],MyImage.shape[1]),None,MyImageCL)
|
|
| 321 |
CLLaunch=BlackHoleCL.EachPixel(queue,(numpy.int32(zImage.shape[0]),numpy.int32(zImage.shape[1])),None,zImageCL,fImageCL)
|
|
| 308 | 322 |
|
| 309 | 323 |
CLLaunch.wait() |
| 310 | 324 |
elapsed = time.time()-start_time |
| 311 | 325 |
print("Elapsed %f: " % elapsed)
|
| 312 | 326 |
|
| 313 |
cl.enqueue_copy(queue,MyImage,MyImageCL).wait() |
|
| 314 |
|
|
| 315 |
print(MyImage) |
|
| 316 |
MyImageCL.release() |
|
| 327 |
cl.enqueue_copy(queue,zImage,zImageCL).wait() |
|
| 328 |
cl.enqueue_copy(queue,fImage,fImageCL).wait() |
|
| 317 | 329 |
|
| 318 |
ImageOutput(MyImage,"TrouNoir") |
|
| 330 |
print(zImage.max()) |
|
| 331 |
print(fImage.max()) |
|
| 332 |
zImageCL.release() |
|
| 333 |
fImageCL.release() |
|
| 334 |
|
|
| 335 |
ImageOutput(zImage,"TrouNoirZ") |
|
| 336 |
ImageOutput(fImage,"TrouNoirF") |
|
| 319 | 337 |
return(elapsed) |
| 320 | 338 |
|
| 321 | 339 |
if __name__=='__main__': |
| ... | ... | |
| 460 | 478 |
# print(Devices,Alu) |
| 461 | 479 |
|
| 462 | 480 |
# MyImage=numpy.where(numpy.random.zeros(Size,Size)>0,1,-1).astype(numpy.float32) |
| 463 |
MyImage=numpy.zeros((Size,Size),dtype=numpy.float32) |
|
| 481 |
zImage=numpy.zeros((Size,Size),dtype=numpy.float32) |
|
| 482 |
fImage=numpy.zeros((Size,Size),dtype=numpy.float32) |
|
| 464 | 483 |
|
| 465 |
print(MyImage)
|
|
| 484 |
print(zImage)
|
|
| 466 | 485 |
|
| 467 |
duration=BlackHole(MyImage,Device)
|
|
| 486 |
duration=BlackHole(zImage,fImage,Device)
|
|
| 468 | 487 |
|
| 469 | 488 |
|
Formats disponibles : Unified diff