From 870236b3c4950762a73247c68023a8dee6e14a7b Mon Sep 17 00:00:00 2001 From: Anthony LaTorre Date: Sun, 12 Jun 2011 21:31:22 -0400 Subject: added some fun models; added some untested code to implement absorption, scattering, reflection, and refraction --- chromaticity.py | 18 ++ ciexyz64_1.csv | 471 +++++++++++++++++++++++++++++++++++++++ geometry.py | 48 ++-- materials.py | 10 +- models/Colbert_HighRes_Brow.STL | Bin 0 -> 13910384 bytes models/Colbert_HighRes_Smile.STL | Bin 0 -> 8699284 bytes models/film.stl | 86 +++++++ models/liberty.stl | Bin 0 -> 853084 bytes photon.py | 29 --- photon_map.py | 86 +++++++ sample.py | 29 +++ solid.py | 19 +- src/intersect.h | 6 + src/kernel.cu | 268 +++++++++++++++++++--- src/materials.h | 45 +++- src/physical_constants.h | 7 + src/random.h | 22 ++ 17 files changed, 1038 insertions(+), 106 deletions(-) create mode 100644 chromaticity.py create mode 100644 ciexyz64_1.csv create mode 100644 models/Colbert_HighRes_Brow.STL create mode 100644 models/Colbert_HighRes_Smile.STL create mode 100644 models/film.stl create mode 100644 models/liberty.stl delete mode 100644 photon.py create mode 100644 photon_map.py create mode 100644 sample.py create mode 100644 src/physical_constants.h create mode 100644 src/random.h diff --git a/chromaticity.py b/chromaticity.py new file mode 100644 index 0000000..18ec28f --- /dev/null +++ b/chromaticity.py @@ -0,0 +1,18 @@ +import numpy as np + +f = open('ciexyz64_1.csv') + +color_map = [] +for line in f: + color_map.append([float(s) for s in line.split(',')]) + +f.close() + +color_map = np.array(color_map) + +def map_wavelength(wavelength): + r = np.interp(wavelength, color_map[:,0], color_map[:,1]) + g = np.interp(wavelength, color_map[:,0], color_map[:,2]) + b = np.interp(wavelength, color_map[:,0], color_map[:,3]) + + return r, g, b diff --git a/ciexyz64_1.csv b/ciexyz64_1.csv new file mode 100644 index 0000000..f437db7 --- /dev/null +++ b/ciexyz64_1.csv @@ -0,0 +1,471 @@ +360,0.000000122200,0.000000013398,0.000000535027 +361,0.000000185138,0.000000020294,0.000000810720 +362,0.000000278830,0.000000030560,0.000001221200 +363,0.000000417470,0.000000045740,0.000001828700 +364,0.000000621330,0.000000068050,0.000002722200 +365,0.000000919270,0.000000100650,0.000004028300 +366,0.000001351980,0.000000147980,0.000005925700 +367,0.000001976540,0.000000216270,0.000008665100 +368,0.000002872500,0.000000314200,0.000012596000 +369,0.000004149500,0.000000453700,0.000018201000 +370,0.000005958600,0.000000651100,0.000026143700 +371,0.000008505600,0.000000928800,0.000037330000 +372,0.000012068600,0.000001317500,0.000052987000 +373,0.000017022600,0.000001857200,0.000074764000 +374,0.000023868000,0.000002602000,0.000104870000 +375,0.000033266000,0.000003625000,0.000146220000 +376,0.000046087000,0.000005019000,0.000202660000 +377,0.000063472000,0.000006907000,0.000279230000 +378,0.000086892000,0.000009449000,0.000382450000 +379,0.000118246000,0.000012848000,0.000520720000 +380,0.000159952000,0.000017364000,0.000704776000 +381,0.000215080000,0.000023327000,0.000948230000 +382,0.000287490000,0.000031150000,0.001268200000 +383,0.000381990000,0.000041350000,0.001686100000 +384,0.000504550000,0.000054560000,0.002228500000 +385,0.000662440000,0.000071560000,0.002927800000 +386,0.000864500000,0.000093300000,0.003823700000 +387,0.001121500000,0.000120870000,0.004964200000 +388,0.001446160000,0.000155640000,0.006406700000 +389,0.001853590000,0.000199200000,0.008219300000 +390,0.002361600000,0.000253400000,0.010482200000 +391,0.002990600000,0.000320200000,0.013289000000 +392,0.003764500000,0.000402400000,0.016747000000 +393,0.004710200000,0.000502300000,0.020980000000 +394,0.005858100000,0.000623200000,0.026127000000 +395,0.007242300000,0.000768500000,0.032344000000 +396,0.008899600000,0.000941700000,0.039802000000 +397,0.010870900000,0.001147800000,0.048691000000 +398,0.013198900000,0.001390300000,0.059210000000 +399,0.015929200000,0.001674000000,0.071576000000 +400,0.019109700000,0.002004400000,0.086010900000 +401,0.022788000000,0.002386000000,0.102740000000 +402,0.027011000000,0.002822000000,0.122000000000 +403,0.031829000000,0.003319000000,0.144020000000 +404,0.037278000000,0.003880000000,0.168990000000 +405,0.043400000000,0.004509000000,0.197120000000 +406,0.050223000000,0.005209000000,0.228570000000 +407,0.057764000000,0.005985000000,0.263470000000 +408,0.066038000000,0.006833000000,0.301900000000 +409,0.075033000000,0.007757000000,0.343870000000 +410,0.084736000000,0.008756000000,0.389366000000 +411,0.095041000000,0.009816000000,0.437970000000 +412,0.105836000000,0.010918000000,0.489220000000 +413,0.117066000000,0.012058000000,0.542900000000 +414,0.128682000000,0.013237000000,0.598810000000 +415,0.140638000000,0.014456000000,0.656760000000 +416,0.152893000000,0.015717000000,0.716580000000 +417,0.165416000000,0.017025000000,0.778120000000 +418,0.178191000000,0.018399000000,0.841310000000 +419,0.191214000000,0.019848000000,0.906110000000 +420,0.204492000000,0.021391000000,0.972542000000 +421,0.217650000000,0.022992000000,1.038900000000 +422,0.230267000000,0.024598000000,1.103100000000 +423,0.242311000000,0.026213000000,1.165100000000 +424,0.253793000000,0.027841000000,1.224900000000 +425,0.264737000000,0.029497000000,1.282500000000 +426,0.275195000000,0.031195000000,1.338200000000 +427,0.285301000000,0.032927000000,1.392600000000 +428,0.295143000000,0.034738000000,1.446100000000 +429,0.304869000000,0.036654000000,1.499400000000 +430,0.314679000000,0.038676000000,1.553480000000 +431,0.324355000000,0.040792000000,1.607200000000 +432,0.333570000000,0.042946000000,1.658900000000 +433,0.342243000000,0.045114000000,1.708200000000 +434,0.350312000000,0.047333000000,1.754800000000 +435,0.357719000000,0.049602000000,1.798500000000 +436,0.364482000000,0.051934000000,1.839200000000 +437,0.370493000000,0.054337000000,1.876600000000 +438,0.375727000000,0.056822000000,1.910500000000 +439,0.380158000000,0.059399000000,1.940800000000 +440,0.383734000000,0.062077000000,1.967280000000 +441,0.386327000000,0.064737000000,1.989100000000 +442,0.387858000000,0.067285000000,2.005700000000 +443,0.388396000000,0.069764000000,2.017400000000 +444,0.387978000000,0.072218000000,2.024400000000 +445,0.386726000000,0.074704000000,2.027300000000 +446,0.384696000000,0.077272000000,2.026400000000 +447,0.382006000000,0.079979000000,2.022300000000 +448,0.378709000000,0.082874000000,2.015300000000 +449,0.374915000000,0.086000000000,2.006000000000 +450,0.370702000000,0.089456000000,1.994800000000 +451,0.366089000000,0.092947000000,1.981400000000 +452,0.361045000000,0.096275000000,1.965300000000 +453,0.355518000000,0.099535000000,1.946400000000 +454,0.349486000000,0.102829000000,1.924800000000 +455,0.342957000000,0.106256000000,1.900700000000 +456,0.335893000000,0.109901000000,1.874100000000 +457,0.328284000000,0.113835000000,1.845100000000 +458,0.320150000000,0.118167000000,1.813900000000 +459,0.311475000000,0.122932000000,1.780600000000 +460,0.302273000000,0.128201000000,1.745370000000 +461,0.292858000000,0.133457000000,1.709100000000 +462,0.283502000000,0.138323000000,1.672300000000 +463,0.274044000000,0.143042000000,1.634700000000 +464,0.264263000000,0.147787000000,1.595600000000 +465,0.254085000000,0.152761000000,1.554900000000 +466,0.243392000000,0.158102000000,1.512200000000 +467,0.232187000000,0.163941000000,1.467300000000 +468,0.220488000000,0.170362000000,1.419900000000 +469,0.208198000000,0.177425000000,1.370000000000 +470,0.195618000000,0.185190000000,1.317560000000 +471,0.183034000000,0.193025000000,1.262400000000 +472,0.170222000000,0.200313000000,1.205000000000 +473,0.157348000000,0.207156000000,1.146600000000 +474,0.144650000000,0.213644000000,1.088000000000 +475,0.132349000000,0.219940000000,1.030200000000 +476,0.120584000000,0.226170000000,0.973830000000 +477,0.109456000000,0.232467000000,0.919430000000 +478,0.099042000000,0.239025000000,0.867460000000 +479,0.089388000000,0.245997000000,0.818280000000 +480,0.080507000000,0.253589000000,0.772125000000 +481,0.072034000000,0.261876000000,0.728290000000 +482,0.063710000000,0.270643000000,0.686040000000 +483,0.055694000000,0.279645000000,0.645530000000 +484,0.048117000000,0.288694000000,0.606850000000 +485,0.041072000000,0.297665000000,0.570060000000 +486,0.034642000000,0.306469000000,0.535220000000 +487,0.028896000000,0.315035000000,0.502340000000 +488,0.023876000000,0.323335000000,0.471400000000 +489,0.019628000000,0.331366000000,0.442390000000 +490,0.016172000000,0.339133000000,0.415254000000 +491,0.013300000000,0.347860000000,0.390024000000 +492,0.010759000000,0.358326000000,0.366399000000 +493,0.008542000000,0.370001000000,0.344015000000 +494,0.006661000000,0.382464000000,0.322689000000 +495,0.005132000000,0.395379000000,0.302356000000 +496,0.003982000000,0.408482000000,0.283036000000 +497,0.003239000000,0.421588000000,0.264816000000 +498,0.002934000000,0.434619000000,0.247848000000 +499,0.003114000000,0.447601000000,0.232318000000 +500,0.003816000000,0.460777000000,0.218502000000 +501,0.005095000000,0.474340000000,0.205851000000 +502,0.006936000000,0.488200000000,0.193596000000 +503,0.009299000000,0.502340000000,0.181736000000 +504,0.012147000000,0.516740000000,0.170281000000 +505,0.015444000000,0.531360000000,0.159249000000 +506,0.019156000000,0.546190000000,0.148673000000 +507,0.023250000000,0.561180000000,0.138609000000 +508,0.027690000000,0.576290000000,0.129096000000 +509,0.032444000000,0.591500000000,0.120215000000 +510,0.037465000000,0.606741000000,0.112044000000 +511,0.042956000000,0.622150000000,0.104710000000 +512,0.049114000000,0.637830000000,0.098196000000 +513,0.055920000000,0.653710000000,0.092361000000 +514,0.063349000000,0.669680000000,0.087088000000 +515,0.071358000000,0.685660000000,0.082248000000 +516,0.079901000000,0.701550000000,0.077744000000 +517,0.088909000000,0.717230000000,0.073456000000 +518,0.098293000000,0.732570000000,0.069268000000 +519,0.107949000000,0.747460000000,0.065060000000 +520,0.117749000000,0.761757000000,0.060709000000 +521,0.127839000000,0.775340000000,0.056457000000 +522,0.138450000000,0.788220000000,0.052609000000 +523,0.149516000000,0.800460000000,0.049122000000 +524,0.161041000000,0.812140000000,0.045954000000 +525,0.172953000000,0.823330000000,0.043050000000 +526,0.185209000000,0.834120000000,0.040368000000 +527,0.197755000000,0.844600000000,0.037839000000 +528,0.210538000000,0.854870000000,0.035384000000 +529,0.223460000000,0.865040000000,0.032949000000 +530,0.236491000000,0.875211000000,0.030451000000 +531,0.249633000000,0.885370000000,0.028029000000 +532,0.262972000000,0.895370000000,0.025862000000 +533,0.276515000000,0.905150000000,0.023920000000 +534,0.290269000000,0.914650000000,0.022174000000 +535,0.304213000000,0.923810000000,0.020584000000 +536,0.318361000000,0.932550000000,0.019127000000 +537,0.332705000000,0.940810000000,0.017740000000 +538,0.347232000000,0.948520000000,0.016403000000 +539,0.361926000000,0.955600000000,0.015064000000 +540,0.376772000000,0.961988000000,0.013676000000 +541,0.391683000000,0.967540000000,0.012308000000 +542,0.406594000000,0.972230000000,0.011056000000 +543,0.421539000000,0.976170000000,0.009915000000 +544,0.436517000000,0.979460000000,0.008872000000 +545,0.451584000000,0.982200000000,0.007918000000 +546,0.466782000000,0.984520000000,0.007030000000 +547,0.482147000000,0.986520000000,0.006223000000 +548,0.497738000000,0.988320000000,0.005453000000 +549,0.513606000000,0.990020000000,0.004714000000 +550,0.529826000000,0.991761000000,0.003988000000 +551,0.546440000000,0.993530000000,0.003289000000 +552,0.563426000000,0.995230000000,0.002646000000 +553,0.580726000000,0.996770000000,0.002063000000 +554,0.598290000000,0.998090000000,0.001533000000 +555,0.616053000000,0.999110000000,0.001091000000 +556,0.633948000000,0.999770000000,0.000711000000 +557,0.651901000000,1.000000000000,0.000407000000 +558,0.669824000000,0.999710000000,0.000184000000 +559,0.687632000000,0.998850000000,0.000047000000 +560,0.705224000000,0.997340000000,0.000000000000 +561,0.722773000000,0.995260000000,0.000000000000 +562,0.740483000000,0.992740000000,0.000000000000 +563,0.758273000000,0.989750000000,0.000000000000 +564,0.776083000000,0.986300000000,0.000000000000 +565,0.793832000000,0.982380000000,0.000000000000 +566,0.811436000000,0.977980000000,0.000000000000 +567,0.828822000000,0.973110000000,0.000000000000 +568,0.845879000000,0.967740000000,0.000000000000 +569,0.862525000000,0.961890000000,0.000000000000 +570,0.878655000000,0.955552000000,0.000000000000 +571,0.894208000000,0.948601000000,0.000000000000 +572,0.909206000000,0.940981000000,0.000000000000 +573,0.923672000000,0.932798000000,0.000000000000 +574,0.937638000000,0.924158000000,0.000000000000 +575,0.951162000000,0.915175000000,0.000000000000 +576,0.964283000000,0.905954000000,0.000000000000 +577,0.977068000000,0.896608000000,0.000000000000 +578,0.989590000000,0.887249000000,0.000000000000 +579,1.001910000000,0.877986000000,0.000000000000 +580,1.014160000000,0.868934000000,0.000000000000 +581,1.026500000000,0.860164000000,0.000000000000 +582,1.038800000000,0.851519000000,0.000000000000 +583,1.051000000000,0.842963000000,0.000000000000 +584,1.062900000000,0.834393000000,0.000000000000 +585,1.074300000000,0.825623000000,0.000000000000 +586,1.085200000000,0.816764000000,0.000000000000 +587,1.095200000000,0.807544000000,0.000000000000 +588,1.104200000000,0.797947000000,0.000000000000 +589,1.112000000000,0.787893000000,0.000000000000 +590,1.118520000000,0.777405000000,0.000000000000 +591,1.123800000000,0.766490000000,0.000000000000 +592,1.128000000000,0.755309000000,0.000000000000 +593,1.131100000000,0.743845000000,0.000000000000 +594,1.133200000000,0.732190000000,0.000000000000 +595,1.134300000000,0.720353000000,0.000000000000 +596,1.134300000000,0.708281000000,0.000000000000 +597,1.133300000000,0.696055000000,0.000000000000 +598,1.131200000000,0.683621000000,0.000000000000 +599,1.128100000000,0.671048000000,0.000000000000 +600,1.123990000000,0.658341000000,0.000000000000 +601,1.118900000000,0.645545000000,0.000000000000 +602,1.112900000000,0.632718000000,0.000000000000 +603,1.105900000000,0.619815000000,0.000000000000 +604,1.098000000000,0.606887000000,0.000000000000 +605,1.089100000000,0.593878000000,0.000000000000 +606,1.079200000000,0.580781000000,0.000000000000 +607,1.068400000000,0.567653000000,0.000000000000 +608,1.056700000000,0.554490000000,0.000000000000 +609,1.044000000000,0.541228000000,0.000000000000 +610,1.030480000000,0.527963000000,0.000000000000 +611,1.016000000000,0.514634000000,0.000000000000 +612,1.000800000000,0.501363000000,0.000000000000 +613,0.984790000000,0.488124000000,0.000000000000 +614,0.968080000000,0.474935000000,0.000000000000 +615,0.950740000000,0.461834000000,0.000000000000 +616,0.932800000000,0.448823000000,0.000000000000 +617,0.914340000000,0.435917000000,0.000000000000 +618,0.895390000000,0.423153000000,0.000000000000 +619,0.876030000000,0.410526000000,0.000000000000 +620,0.856297000000,0.398057000000,0.000000000000 +621,0.836350000000,0.385835000000,0.000000000000 +622,0.816290000000,0.373951000000,0.000000000000 +623,0.796050000000,0.362311000000,0.000000000000 +624,0.775610000000,0.350863000000,0.000000000000 +625,0.754930000000,0.339554000000,0.000000000000 +626,0.733990000000,0.328309000000,0.000000000000 +627,0.712780000000,0.317118000000,0.000000000000 +628,0.691290000000,0.305936000000,0.000000000000 +629,0.669520000000,0.294737000000,0.000000000000 +630,0.647467000000,0.283493000000,0.000000000000 +631,0.625110000000,0.272222000000,0.000000000000 +632,0.602520000000,0.260990000000,0.000000000000 +633,0.579890000000,0.249877000000,0.000000000000 +634,0.557370000000,0.238946000000,0.000000000000 +635,0.535110000000,0.228254000000,0.000000000000 +636,0.513240000000,0.217853000000,0.000000000000 +637,0.491860000000,0.207780000000,0.000000000000 +638,0.471080000000,0.198072000000,0.000000000000 +639,0.450960000000,0.188748000000,0.000000000000 +640,0.431567000000,0.179828000000,0.000000000000 +641,0.412870000000,0.171285000000,0.000000000000 +642,0.394750000000,0.163059000000,0.000000000000 +643,0.377210000000,0.155151000000,0.000000000000 +644,0.360190000000,0.147535000000,0.000000000000 +645,0.343690000000,0.140211000000,0.000000000000 +646,0.327690000000,0.133170000000,0.000000000000 +647,0.312170000000,0.126400000000,0.000000000000 +648,0.297110000000,0.119892000000,0.000000000000 +649,0.282500000000,0.113640000000,0.000000000000 +650,0.268329000000,0.107633000000,0.000000000000 +651,0.254590000000,0.101870000000,0.000000000000 +652,0.241300000000,0.096347000000,0.000000000000 +653,0.228480000000,0.091063000000,0.000000000000 +654,0.216140000000,0.086010000000,0.000000000000 +655,0.204300000000,0.081187000000,0.000000000000 +656,0.192950000000,0.076583000000,0.000000000000 +657,0.182110000000,0.072198000000,0.000000000000 +658,0.171770000000,0.068024000000,0.000000000000 +659,0.161920000000,0.064052000000,0.000000000000 +660,0.152568000000,0.060281000000,0.000000000000 +661,0.143670000000,0.056697000000,0.000000000000 +662,0.135200000000,0.053292000000,0.000000000000 +663,0.127130000000,0.050059000000,0.000000000000 +664,0.119480000000,0.046998000000,0.000000000000 +665,0.112210000000,0.044096000000,0.000000000000 +666,0.105310000000,0.041345000000,0.000000000000 +667,0.098786000000,0.038750700000,0.000000000000 +668,0.092610000000,0.036297800000,0.000000000000 +669,0.086773000000,0.033983200000,0.000000000000 +670,0.081260600000,0.031800400000,0.000000000000 +671,0.076048000000,0.029739500000,0.000000000000 +672,0.071114000000,0.027791800000,0.000000000000 +673,0.066454000000,0.025955100000,0.000000000000 +674,0.062062000000,0.024226300000,0.000000000000 +675,0.057930000000,0.022601700000,0.000000000000 +676,0.054050000000,0.021077900000,0.000000000000 +677,0.050412000000,0.019650500000,0.000000000000 +678,0.047006000000,0.018315300000,0.000000000000 +679,0.043823000000,0.017068600000,0.000000000000 +680,0.040850800000,0.015905100000,0.000000000000 +681,0.038072000000,0.014818300000,0.000000000000 +682,0.035468000000,0.013800800000,0.000000000000 +683,0.033031000000,0.012849500000,0.000000000000 +684,0.030753000000,0.011960700000,0.000000000000 +685,0.028623000000,0.011130300000,0.000000000000 +686,0.026635000000,0.010355500000,0.000000000000 +687,0.024781000000,0.009633200000,0.000000000000 +688,0.023052000000,0.008959900000,0.000000000000 +689,0.021441000000,0.008332400000,0.000000000000 +690,0.019941300000,0.007748800000,0.000000000000 +691,0.018544000000,0.007204600000,0.000000000000 +692,0.017241000000,0.006697500000,0.000000000000 +693,0.016027000000,0.006225100000,0.000000000000 +694,0.014896000000,0.005785000000,0.000000000000 +695,0.013842000000,0.005375100000,0.000000000000 +696,0.012862000000,0.004994100000,0.000000000000 +697,0.011949000000,0.004639200000,0.000000000000 +698,0.011100000000,0.004309300000,0.000000000000 +699,0.010311000000,0.004002800000,0.000000000000 +700,0.009576880000,0.003717740000,0.000000000000 +701,0.008894000000,0.003452620000,0.000000000000 +702,0.008258100000,0.003205830000,0.000000000000 +703,0.007666400000,0.002976230000,0.000000000000 +704,0.007116300000,0.002762810000,0.000000000000 +705,0.006605200000,0.002564560000,0.000000000000 +706,0.006130600000,0.002380480000,0.000000000000 +707,0.005690300000,0.002209710000,0.000000000000 +708,0.005281900000,0.002051320000,0.000000000000 +709,0.004903300000,0.001904490000,0.000000000000 +710,0.004552630000,0.001768470000,0.000000000000 +711,0.004227500000,0.001642360000,0.000000000000 +712,0.003925800000,0.001525350000,0.000000000000 +713,0.003645700000,0.001416720000,0.000000000000 +714,0.003385900000,0.001315950000,0.000000000000 +715,0.003144700000,0.001222390000,0.000000000000 +716,0.002920800000,0.001135550000,0.000000000000 +717,0.002713000000,0.001054940000,0.000000000000 +718,0.002520200000,0.000980140000,0.000000000000 +719,0.002341100000,0.000910660000,0.000000000000 +720,0.002174960000,0.000846190000,0.000000000000 +721,0.002020600000,0.000786290000,0.000000000000 +722,0.001877300000,0.000730680000,0.000000000000 +723,0.001744100000,0.000678990000,0.000000000000 +724,0.001620500000,0.000631010000,0.000000000000 +725,0.001505700000,0.000586440000,0.000000000000 +726,0.001399200000,0.000545110000,0.000000000000 +727,0.001300400000,0.000506720000,0.000000000000 +728,0.001208700000,0.000471110000,0.000000000000 +729,0.001123600000,0.000438050000,0.000000000000 +730,0.001044760000,0.000407410000,0.000000000000 +731,0.000971560000,0.000378962000,0.000000000000 +732,0.000903600000,0.000352543000,0.000000000000 +733,0.000840480000,0.000328001000,0.000000000000 +734,0.000781870000,0.000305208000,0.000000000000 +735,0.000727450000,0.000284041000,0.000000000000 +736,0.000676900000,0.000264375000,0.000000000000 +737,0.000629960000,0.000246109000,0.000000000000 +738,0.000586370000,0.000229143000,0.000000000000 +739,0.000545870000,0.000213376000,0.000000000000 +740,0.000508258000,0.000198730000,0.000000000000 +741,0.000473300000,0.000185115000,0.000000000000 +742,0.000440800000,0.000172454000,0.000000000000 +743,0.000410580000,0.000160678000,0.000000000000 +744,0.000382490000,0.000149730000,0.000000000000 +745,0.000356380000,0.000139550000,0.000000000000 +746,0.000332110000,0.000130086000,0.000000000000 +747,0.000309550000,0.000121290000,0.000000000000 +748,0.000288580000,0.000113106000,0.000000000000 +749,0.000269090000,0.000105501000,0.000000000000 +750,0.000250969000,0.000098428000,0.000000000000 +751,0.000234130000,0.000091853000,0.000000000000 +752,0.000218470000,0.000085738000,0.000000000000 +753,0.000203910000,0.000080048000,0.000000000000 +754,0.000190350000,0.000074751000,0.000000000000 +755,0.000177730000,0.000069819000,0.000000000000 +756,0.000165970000,0.000065222000,0.000000000000 +757,0.000155020000,0.000060939000,0.000000000000 +758,0.000144800000,0.000056942000,0.000000000000 +759,0.000135280000,0.000053217000,0.000000000000 +760,0.000126390000,0.000049737000,0.000000000000 +761,0.000118100000,0.000046491000,0.000000000000 +762,0.000110370000,0.000043464000,0.000000000000 +763,0.000103150000,0.000040635000,0.000000000000 +764,0.000096427000,0.000038000000,0.000000000000 +765,0.000090151000,0.000035540500,0.000000000000 +766,0.000084294000,0.000033244800,0.000000000000 +767,0.000078830000,0.000031100600,0.000000000000 +768,0.000073729000,0.000029099000,0.000000000000 +769,0.000068969000,0.000027230700,0.000000000000 +770,0.000064525800,0.000025486000,0.000000000000 +771,0.000060376000,0.000023856100,0.000000000000 +772,0.000056500000,0.000022333200,0.000000000000 +773,0.000052880000,0.000020910400,0.000000000000 +774,0.000049498000,0.000019580800,0.000000000000 +775,0.000046339000,0.000018338400,0.000000000000 +776,0.000043389000,0.000017177700,0.000000000000 +777,0.000040634000,0.000016093400,0.000000000000 +778,0.000038060000,0.000015080000,0.000000000000 +779,0.000035657000,0.000014133600,0.000000000000 +780,0.000033411700,0.000013249000,0.000000000000 +781,0.000031315000,0.000012422600,0.000000000000 +782,0.000029355000,0.000011649900,0.000000000000 +783,0.000027524000,0.000010927700,0.000000000000 +784,0.000025811000,0.000010251900,0.000000000000 +785,0.000024209000,0.000009619600,0.000000000000 +786,0.000022711000,0.000009028100,0.000000000000 +787,0.000021308000,0.000008474000,0.000000000000 +788,0.000019994000,0.000007954800,0.000000000000 +789,0.000018764000,0.000007468600,0.000000000000 +790,0.000017611500,0.000007012800,0.000000000000 +791,0.000016532000,0.000006585800,0.000000000000 +792,0.000015521000,0.000006185700,0.000000000000 +793,0.000014574000,0.000005810700,0.000000000000 +794,0.000013686000,0.000005459000,0.000000000000 +795,0.000012855000,0.000005129800,0.000000000000 +796,0.000012075000,0.000004820600,0.000000000000 +797,0.000011345000,0.000004531200,0.000000000000 +798,0.000010659000,0.000004259100,0.000000000000 +799,0.000010017000,0.000004004200,0.000000000000 +800,0.000009413630,0.000003764730,0.000000000000 +801,0.000008847900,0.000003539950,0.000000000000 +802,0.000008317100,0.000003329140,0.000000000000 +803,0.000007819000,0.000003131150,0.000000000000 +804,0.000007351600,0.000002945290,0.000000000000 +805,0.000006913000,0.000002770810,0.000000000000 +806,0.000006501500,0.000002607050,0.000000000000 +807,0.000006115300,0.000002453290,0.000000000000 +808,0.000005752900,0.000002308940,0.000000000000 +809,0.000005412700,0.000002173380,0.000000000000 +810,0.000005093470,0.000002046130,0.000000000000 +811,0.000004793800,0.000001926620,0.000000000000 +812,0.000004512500,0.000001814400,0.000000000000 +813,0.000004248300,0.000001708950,0.000000000000 +814,0.000004000200,0.000001609880,0.000000000000 +815,0.000003767100,0.000001516770,0.000000000000 +816,0.000003548000,0.000001429210,0.000000000000 +817,0.000003342100,0.000001346860,0.000000000000 +818,0.000003148500,0.000001269450,0.000000000000 +819,0.000002966500,0.000001196620,0.000000000000 +820,0.000002795310,0.000001128090,0.000000000000 +821,0.000002634500,0.000001063680,0.000000000000 +822,0.000002483400,0.000001003130,0.000000000000 +823,0.000002341400,0.000000946220,0.000000000000 +824,0.000002207800,0.000000892630,0.000000000000 +825,0.000002082000,0.000000842160,0.000000000000 +826,0.000001963600,0.000000794640,0.000000000000 +827,0.000001851900,0.000000749780,0.000000000000 +828,0.000001746500,0.000000707440,0.000000000000 +829,0.000001647100,0.000000667480,0.000000000000 +830,0.000001553140,0.000000629700,0.000000000000 \ No newline at end of file diff --git a/geometry.py b/geometry.py index 8af8d49..8dee8c9 100644 --- a/geometry.py +++ b/geometry.py @@ -106,25 +106,17 @@ class Geometry(object): else: self.material2_index[i] = -1 - surface1 = np.concatenate([solid.surface1 for solid in self.solids]) - surface2 = np.concatenate([solid.surface2 for solid in self.solids]) + surface = np.concatenate([solid.surface for solid in self.solids]) - self.surfaces = list(np.unique(np.concatenate([surface1, surface2]))) + self.surfaces = list(np.unique(surface)) - self.surface1_index = np.empty(len(self.mesh), dtype=np.int32) - self.surface2_index = np.empty(len(self.mesh), dtype=np.int32) + self.surface_index = np.empty(len(self.mesh), dtype=np.int32) - for i, surface in enumerate(surface1[reorder]): + for i, surface in enumerate(surface[reorder]): if surface is not None: - self.surface1_index[i] = self.surfaces.index(surface) + self.surface_index[i] = self.surfaces.index(surface) else: - self.surface1_index[i] = -1 - - for i, surface in enumerate(surface2[reorder]): - if surface is not None: - self.surface2_index[i] = self.surfaces.index(surface) - else: - self.surface2_index[i] = -1 + self.surface_index[i] = -1 self.colors = np.concatenate([solid.color for solid in self.solids]) @@ -210,6 +202,10 @@ class Geometry(object): of the texture references. """ + set_wavelength_range = module.get_function('set_wavelength_range') + + set_wavelength_range(wavelengths[0], wavelengths[-1], wavelengths[1]-wavelengths[0], block=(1,1,1), grid=(1,1)) + set_material = module.get_function('set_material') for i, material in enumerate(self.materials): if material is None: @@ -223,7 +219,7 @@ class Geometry(object): material.absorption_length_gpu = cuda.to_device(absorption_length) material.scattering_length_gpu = cuda.to_device(scattering_length) - set_material(i, material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1)) + set_material(np.int32(i), material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1)) set_surface = module.get_function('set_surface') for i, surface in enumerate(self.surfaces): @@ -232,33 +228,29 @@ class Geometry(object): absorption = np.interp(wavelengths, surface.absorption[:,0], surface.absorption[:,1]).astype(np.float32) reflection_diffuse = np.interp(wavelengths, surface.reflection_diffuse[:,0], surface.reflection_diffuse[:,1]).astype(np.float32) - reflection_specular = np.interp(wavelengths, surface.refelection_specular[:,0], surface.reflection_specular[:,1]).astype(np.float32) + reflection_specular = np.interp(wavelengths, surface.reflection_specular[:,0], surface.reflection_specular[:,1]).astype(np.float32) surface.absorption_gpu = cuda.to_device(absorption) surface.reflection_diffuse_gpu = cuda.to_device(reflection_diffuse) surface.reflection_specular_gpu = cuda.to_device(reflection_specular) - set_surface(i, surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu) + set_surface(np.int32(i), surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu, block=(1,1,1), grid=(1,1)) - material1_index_tex = module.get_texref('material1_index') - material2_index_tex = module.get_texref('material2_index') - surface1_index_tex = module.get_texref('surface1_index') - surface2_index_tex = module.get_texref('surface2_index') + material1_index_tex = module.get_texref('material1_lookup') + material2_index_tex = module.get_texref('material2_lookup') + surface_index_tex = module.get_texref('surface_lookup') material1_index_gpu = cuda.to_device(self.material1_index) material2_index_gpu = cuda.to_device(self.material2_index) - surface1_index_gpu = cuda.to_device(self.surface1_index) - surface2_index_gpu = cuda.to_device(self.surface2_index) + surface_index_gpu = cuda.to_device(self.surface_index) material1_index_tex.set_address(material1_index_gpu, self.material1_index.nbytes) material2_index_tex.set_address(material2_index_gpu, self.material2_index.nbytes) - surface1_index_tex.set_address(surface1_index_gpu, self.surface1_index.nbytes) - surface2_index_tex.set_address(surface2_index_gpu, self.surface2_index.nbytes) + surface_index_tex.set_address(surface_index_gpu, self.surface_index.nbytes) material1_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) material2_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - surface1_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - surface2_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) + surface_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) vertices = np.empty(len(self.mesh.vertices), dtype=gpuarray.vec.float4) vertices['x'] = self.mesh.vertices[:,0] @@ -274,7 +266,7 @@ class Geometry(object): if color: triangles['w'] = self.colors else: - triangles['w'] = (self.material1_index << 24) | (self.material2_index << 16) | (self.surface1_index << 8) | self.surface2_index + triangles['w'] = (self.material1_index << 24) | (self.material2_index << 16) | (self.surface_index << 8) lower_bounds = np.empty(self.lower_bounds.shape[0], dtype=gpuarray.vec.float4) lower_bounds['x'] = self.lower_bounds[:,0] diff --git a/materials.py b/materials.py index dc28a8b..273a4bb 100644 --- a/materials.py +++ b/materials.py @@ -1,6 +1,6 @@ import numpy as np -wavelengths = np.arange(200, 810, 20) +wavelengths = np.arange(200, 810, 20).astype(np.float32) class Material(object): """Material optical properties.""" @@ -68,10 +68,10 @@ r1408_photocathode = Surface('photocathode') r1408_photocathode.absorption = \ np.array([(60.0, 0.0), (200.0, 0.0), (230.0, 0.0), (231.0, 0.0), (232.0, 0.0), (233.0, 0.0), (234.0, 0.0), (235.0, 0.0), (236.0, 0.0), (237.0, 0.0), (238.0, 0.0), (239.0, 0.0), (240.0, 0.0), (241.0, 0.00049089999999999995), (242.0, 0.0009817999999999999), (243.0, 0.0014729999999999999), (244.0, 0.001964), (245.0, 0.002454), (246.0, 0.0029450000000000001), (247.0, 0.0034359999999999998), (248.0, 0.0039269999999999999), (249.0, 0.0044180000000000001), (250.0, 0.0049090000000000002), (251.0, 0.0055100000000000001), (252.0, 0.0061110000000000001), (253.0, 0.0067130000000000002), (254.0, 0.0073140000000000002), (255.0, 0.0079150000000000002), (256.0, 0.0085159999999999993), (257.0, 0.0091170000000000001), (258.0, 0.0097190000000000002), (259.0, 0.010319999999999999), (260.0, 0.010919999999999999), (261.0, 0.01269), (262.0, 0.014460000000000001), (263.0, 0.016219999999999998), (264.0, 0.017989999999999999), (265.0, 0.01976), (266.0, 0.021530000000000001), (267.0, 0.023290000000000002), (268.0, 0.025059999999999999), (269.0, 0.02683), (270.0, 0.0286), (271.0, 0.031690000000000003), (272.0, 0.034779999999999998), (273.0, 0.037879999999999997), (274.0, 0.040969999999999999), (275.0, 0.044060000000000002), (276.0, 0.047160000000000001), (277.0, 0.050250000000000003), (278.0, 0.053339999999999999), (279.0, 0.056439999999999997), (280.0, 0.05953), (281.0, 0.062950000000000006), (282.0, 0.066360000000000002), (283.0, 0.069779999999999995), (284.0, 0.073200000000000001), (285.0, 0.076609999999999998), (286.0, 0.080030000000000004), (287.0, 0.08344), (288.0, 0.086860000000000007), (289.0, 0.090279999999999999), (290.0, 0.093689999999999996), (291.0, 0.096799999999999997), (292.0, 0.099909999999999999), (293.0, 0.10299999999999999), (294.0, 0.1061), (295.0, 0.10920000000000001), (296.0, 0.1123), (297.0, 0.1154), (298.0, 0.11849999999999999), (299.0, 0.1217), (300.0, 0.12479999999999999), (301.0, 0.1275), (302.0, 0.13020000000000001), (303.0, 0.13300000000000001), (304.0, 0.13569999999999999), (305.0, 0.13850000000000001), (306.0, 0.14119999999999999), (307.0, 0.1439), (308.0, 0.1467), (309.0, 0.14940000000000001), (310.0, 0.1522), (311.0, 0.15409999999999999), (312.0, 0.15609999999999999), (313.0, 0.158), (314.0, 0.16), (315.0, 0.16189999999999999), (316.0, 0.16389999999999999), (317.0, 0.1658), (318.0, 0.1678), (319.0, 0.16969999999999999), (320.0, 0.17169999999999999), (321.0, 0.17249999999999999), (322.0, 0.1734), (323.0, 0.17419999999999999), (324.0, 0.17510000000000001), (325.0, 0.1759), (326.0, 0.1767), (327.0, 0.17760000000000001), (328.0, 0.1784), (329.0, 0.17929999999999999), (330.0, 0.18010000000000001), (331.0, 0.18099999999999999), (332.0, 0.18190000000000001), (333.0, 0.1827), (334.0, 0.18360000000000001), (335.0, 0.18440000000000001), (336.0, 0.18529999999999999), (337.0, 0.1862), (338.0, 0.187), (339.0, 0.18790000000000001), (340.0, 0.1888), (341.0, 0.189), (342.0, 0.1893), (343.0, 0.18959999999999999), (344.0, 0.1898), (345.0, 0.19009999999999999), (346.0, 0.1903), (347.0, 0.19059999999999999), (348.0, 0.19089999999999999), (349.0, 0.19109999999999999), (350.0, 0.19139999999999999), (351.0, 0.19170000000000001), (352.0, 0.19209999999999999), (353.0, 0.19239999999999999), (354.0, 0.1928), (355.0, 0.19309999999999999), (356.0, 0.19350000000000001), (357.0, 0.1938), (358.0, 0.19420000000000001), (359.0, 0.19450000000000001), (360.0, 0.19489999999999999), (361.0, 0.19520000000000001), (362.0, 0.19539999999999999), (363.0, 0.19570000000000001), (364.0, 0.19600000000000001), (365.0, 0.1963), (366.0, 0.1966), (367.0, 0.19689999999999999), (368.0, 0.19719999999999999), (369.0, 0.19750000000000001), (370.0, 0.1978), (371.0, 0.1978), (372.0, 0.1978), (373.0, 0.1978), (374.0, 0.19789999999999999), (375.0, 0.19789999999999999), (376.0, 0.19789999999999999), (377.0, 0.19789999999999999), (378.0, 0.19800000000000001), (379.0, 0.19800000000000001), (380.0, 0.19800000000000001), (381.0, 0.1983), (382.0, 0.19850000000000001), (383.0, 0.1988), (384.0, 0.19900000000000001), (385.0, 0.1993), (386.0, 0.1996), (387.0, 0.19980000000000001), (388.0, 0.2001), (389.0, 0.20030000000000001), (390.0, 0.2006), (391.0, 0.20080000000000001), (392.0, 0.20100000000000001), (393.0, 0.20130000000000001), (394.0, 0.20150000000000001), (395.0, 0.20169999999999999), (396.0, 0.2019), (397.0, 0.20219999999999999), (398.0, 0.2024), (399.0, 0.2026), (400.0, 0.2029), (401.0, 0.20330000000000001), (402.0, 0.20380000000000001), (403.0, 0.20419999999999999), (404.0, 0.20469999999999999), (405.0, 0.20519999999999999), (406.0, 0.2056), (407.0, 0.20610000000000001), (408.0, 0.20649999999999999), (409.0, 0.20699999999999999), (410.0, 0.20749999999999999), (411.0, 0.20780000000000001), (412.0, 0.2082), (413.0, 0.20849999999999999), (414.0, 0.2089), (415.0, 0.2092), (416.0, 0.20960000000000001), (417.0, 0.2099), (418.0, 0.21029999999999999), (419.0, 0.21060000000000001), (420.0, 0.21099999999999999), (421.0, 0.21129999999999999), (422.0, 0.21160000000000001), (423.0, 0.21190000000000001), (424.0, 0.21229999999999999), (425.0, 0.21260000000000001), (426.0, 0.21290000000000001), (427.0, 0.2132), (428.0, 0.2135), (429.0, 0.21390000000000001), (430.0, 0.2142), (431.0, 0.2142), (432.0, 0.21429999999999999), (433.0, 0.21440000000000001), (434.0, 0.21440000000000001), (435.0, 0.2145), (436.0, 0.21460000000000001), (437.0, 0.21460000000000001), (438.0, 0.2147), (439.0, 0.21479999999999999), (440.0, 0.21479999999999999), (441.0, 0.2147), (442.0, 0.2145), (443.0, 0.21440000000000001), (444.0, 0.2142), (445.0, 0.21410000000000001), (446.0, 0.21390000000000001), (447.0, 0.21379999999999999), (448.0, 0.21360000000000001), (449.0, 0.2135), (450.0, 0.21329999999999999), (451.0, 0.21310000000000001), (452.0, 0.21279999999999999), (453.0, 0.21260000000000001), (454.0, 0.21240000000000001), (455.0, 0.21210000000000001), (456.0, 0.21190000000000001), (457.0, 0.2117), (458.0, 0.2114), (459.0, 0.2112), (460.0, 0.21099999999999999), (461.0, 0.21029999999999999), (462.0, 0.2097), (463.0, 0.20899999999999999), (464.0, 0.2084), (465.0, 0.2077), (466.0, 0.20710000000000001), (467.0, 0.2064), (468.0, 0.20580000000000001), (469.0, 0.2051), (470.0, 0.20449999999999999), (471.0, 0.20380000000000001), (472.0, 0.20319999999999999), (473.0, 0.20250000000000001), (474.0, 0.20180000000000001), (475.0, 0.20119999999999999), (476.0, 0.20050000000000001), (477.0, 0.19989999999999999), (478.0, 0.19919999999999999), (479.0, 0.1986), (480.0, 0.19789999999999999), (481.0, 0.19750000000000001), (482.0, 0.1971), (483.0, 0.19670000000000001), (484.0, 0.19620000000000001), (485.0, 0.1958), (486.0, 0.19539999999999999), (487.0, 0.19500000000000001), (488.0, 0.1946), (489.0, 0.19409999999999999), (490.0, 0.19370000000000001), (491.0, 0.1933), (492.0, 0.19289999999999999), (493.0, 0.1925), (494.0, 0.19209999999999999), (495.0, 0.19170000000000001), (496.0, 0.1913), (497.0, 0.1908), (498.0, 0.19040000000000001), (499.0, 0.19), (500.0, 0.18959999999999999), (501.0, 0.18920000000000001), (502.0, 0.18870000000000001), (503.0, 0.1883), (504.0, 0.18790000000000001), (505.0, 0.18740000000000001), (506.0, 0.187), (507.0, 0.18659999999999999), (508.0, 0.18609999999999999), (509.0, 0.1857), (510.0, 0.1852), (511.0, 0.1835), (512.0, 0.18179999999999999), (513.0, 0.18010000000000001), (514.0, 0.1784), (515.0, 0.17660000000000001), (516.0, 0.1749), (517.0, 0.17319999999999999), (518.0, 0.17150000000000001), (519.0, 0.16980000000000001), (520.0, 0.1681), (521.0, 0.1651), (522.0, 0.16220000000000001), (523.0, 0.1593), (524.0, 0.15640000000000001), (525.0, 0.1535), (526.0, 0.15060000000000001), (527.0, 0.1477), (528.0, 0.14480000000000001), (529.0, 0.1419), (530.0, 0.13900000000000001), (531.0, 0.1366), (532.0, 0.13420000000000001), (533.0, 0.1318), (534.0, 0.12939999999999999), (535.0, 0.127), (536.0, 0.1246), (537.0, 0.1222), (538.0, 0.1198), (539.0, 0.1174), (540.0, 0.115), (541.0, 0.1143), (542.0, 0.11360000000000001), (543.0, 0.1129), (544.0, 0.11219999999999999), (545.0, 0.1115), (546.0, 0.11070000000000001), (547.0, 0.11), (548.0, 0.10929999999999999), (549.0, 0.1086), (550.0, 0.1079), (551.0, 0.1062), (552.0, 0.1046), (553.0, 0.10290000000000001), (554.0, 0.1013), (555.0, 0.099650000000000002), (556.0, 0.09801), (557.0, 0.096360000000000001), (558.0, 0.094719999999999999), (559.0, 0.093079999999999996), (560.0, 0.091429999999999997), (561.0, 0.0906), (562.0, 0.089770000000000003), (563.0, 0.088950000000000001), (564.0, 0.088120000000000004), (565.0, 0.087290000000000006), (566.0, 0.086459999999999995), (567.0, 0.085629999999999998), (568.0, 0.0848), (569.0, 0.083979999999999999), (570.0, 0.083150000000000002), (571.0, 0.082239999999999994), (572.0, 0.08133), (573.0, 0.080409999999999995), (574.0, 0.079500000000000001), (575.0, 0.078589999999999993), (576.0, 0.077679999999999999), (577.0, 0.076770000000000005), (578.0, 0.075859999999999997), (579.0, 0.074950000000000003), (580.0, 0.074039999999999995), (581.0, 0.072969999999999993), (582.0, 0.071900000000000006), (583.0, 0.070830000000000004), (584.0, 0.069760000000000003), (585.0, 0.068690000000000001), (586.0, 0.06762), (587.0, 0.066549999999999998), (588.0, 0.065479999999999997), (589.0, 0.064409999999999995), (590.0, 0.063339999999999994), (591.0, 0.062230000000000001), (592.0, 0.061109999999999998), (593.0, 0.059990000000000002), (594.0, 0.058869999999999999), (595.0, 0.057750000000000003), (596.0, 0.05663), (597.0, 0.05552), (598.0, 0.054399999999999997), (599.0, 0.053280000000000001), (600.0, 0.052159999999999998), (601.0, 0.05092), (602.0, 0.049669999999999999), (603.0, 0.048430000000000001), (604.0, 0.047190000000000003), (605.0, 0.045940000000000002), (606.0, 0.044699999999999997), (607.0, 0.043450000000000003), (608.0, 0.042209999999999998), (609.0, 0.040969999999999999), (610.0, 0.039719999999999998), (611.0, 0.03884), (612.0, 0.037949999999999998), (613.0, 0.037069999999999999), (614.0, 0.036179999999999997), (615.0, 0.035299999999999998), (616.0, 0.034410000000000003), (617.0, 0.033529999999999997), (618.0, 0.032640000000000002), (619.0, 0.031759999999999997), (620.0, 0.030870000000000002), (621.0, 0.029960000000000001), (622.0, 0.029059999999999999), (623.0, 0.028150000000000001), (624.0, 0.02724), (625.0, 0.026339999999999999), (626.0, 0.025430000000000001), (627.0, 0.02452), (628.0, 0.023609999999999999), (629.0, 0.022710000000000001), (630.0, 0.0218), (631.0, 0.021090000000000001), (632.0, 0.020379999999999999), (633.0, 0.01967), (634.0, 0.018970000000000001), (635.0, 0.018259999999999998), (636.0, 0.01755), (637.0, 0.016840000000000001), (638.0, 0.016129999999999999), (639.0, 0.01542), (640.0, 0.014710000000000001), (641.0, 0.014120000000000001), (642.0, 0.01353), (643.0, 0.012930000000000001), (644.0, 0.01234), (645.0, 0.01175), (646.0, 0.01115), (647.0, 0.01056), (648.0, 0.0099670000000000002), (649.0, 0.0093740000000000004), (650.0, 0.0087810000000000006), (651.0, 0.0084159999999999999), (652.0, 0.0080520000000000001), (653.0, 0.0076870000000000003), (654.0, 0.0073229999999999996), (655.0, 0.0069579999999999998), (656.0, 0.006594), (657.0, 0.0062290000000000002), (658.0, 0.0058650000000000004), (659.0, 0.0054999999999999997), (660.0, 0.0051359999999999999), (661.0, 0.0048989999999999997), (662.0, 0.0046629999999999996), (663.0, 0.0044270000000000004), (664.0, 0.0041910000000000003), (665.0, 0.0039550000000000002), (666.0, 0.003718), (667.0, 0.0034819999999999999), (668.0, 0.0032460000000000002), (669.0, 0.0030100000000000001), (670.0, 0.002774), (671.0, 0.00265), (672.0, 0.0025270000000000002), (673.0, 0.0024039999999999999), (674.0, 0.002281), (675.0, 0.0021570000000000001), (676.0, 0.0020339999999999998), (677.0, 0.0019109999999999999), (678.0, 0.0017880000000000001), (679.0, 0.001665), (680.0, 0.0015410000000000001), (681.0, 0.001459), (682.0, 0.001377), (683.0, 0.0012949999999999999), (684.0, 0.0012130000000000001), (685.0, 0.0011310000000000001), (686.0, 0.001049), (687.0, 0.00096650000000000002), (688.0, 0.00088440000000000003), (689.0, 0.00080230000000000004), (690.0, 0.00072009999999999999), (691.0, 0.00069439999999999997), (692.0, 0.00066859999999999999), (693.0, 0.00064289999999999996), (694.0, 0.00061709999999999998), (695.0, 0.00059139999999999996), (696.0, 0.00056559999999999998), (697.0, 0.00053989999999999995), (698.0, 0.00051409999999999997), (699.0, 0.00048840000000000005), (700.0, 0.00046260000000000002), (800.0, 0.0)]) -r1408_photocathode.specular_reflection = \ +r1408_photocathode.reflection_specular = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) -r1408_photocathode.diffuse_reflection = \ +r1408_photocathode.reflection_diffuse = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) aluminum = Surface('aluminum') @@ -79,8 +79,8 @@ aluminum = Surface('aluminum') aluminum.absorption = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) -aluminum.specular_reflection = \ +aluminum.reflection_specular = \ np.array(zip(wavelengths, np.ones(len(wavelengths)))) -aluminum.diffuse_reflection = \ +aluminum.reflection_diffuse = \ np.array(zip(wavelengths, np.zeros(len(wavelengths)))) diff --git a/models/Colbert_HighRes_Brow.STL b/models/Colbert_HighRes_Brow.STL new file mode 100644 index 0000000..784c877 Binary files /dev/null and b/models/Colbert_HighRes_Brow.STL differ diff --git a/models/Colbert_HighRes_Smile.STL b/models/Colbert_HighRes_Smile.STL new file mode 100644 index 0000000..9cf1864 Binary files /dev/null and b/models/Colbert_HighRes_Smile.STL differ diff --git a/models/film.stl b/models/film.stl new file mode 100644 index 0000000..30c0f65 --- /dev/null +++ b/models/film.stl @@ -0,0 +1,86 @@ +solid OpenSCAD_Model + facet normal -1.000000 0.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal -1.000000 -0.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 0.000000 1.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex 0.005000 -1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 -0.000000 1.000000 + outer loop + vertex -0.005000 1.000000 1.000000 + vertex -0.005000 -1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal -0.000000 -1.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 -1.000000 0.000000 + outer loop + vertex -0.005000 -1.000000 1.000000 + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 -1.000000 1.000000 + endloop + endfacet + facet normal 0.000000 0.000000 -1.000000 + outer loop + vertex -0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex -0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal -0.000000 0.000000 -1.000000 + outer loop + vertex -0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex 0.005000 -1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 1.000000 -0.000000 + outer loop + vertex -0.005000 1.000000 1.000000 + vertex 0.005000 1.000000 1.000000 + vertex -0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 0.000000 1.000000 0.000000 + outer loop + vertex -0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + vertex 0.005000 1.000000 -1.000000 + endloop + endfacet + facet normal 1.000000 0.000000 0.000000 + outer loop + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet + facet normal 1.000000 0.000000 -0.000000 + outer loop + vertex 0.005000 -1.000000 1.000000 + vertex 0.005000 -1.000000 -1.000000 + vertex 0.005000 1.000000 1.000000 + endloop + endfacet +endsolid OpenSCAD_Model diff --git a/models/liberty.stl b/models/liberty.stl new file mode 100644 index 0000000..89c7e87 Binary files /dev/null and b/models/liberty.stl differ diff --git a/photon.py b/photon.py deleted file mode 100644 index b494384..0000000 --- a/photon.py +++ /dev/null @@ -1,29 +0,0 @@ -import numpy as np - -def uniform_sphere(size=None, dtype=np.double): - """ - Generate random points isotropically distributed across the unit sphere. - - Args: - - size: int, *optional* - Number of points to generate. If no size is specified, a single - point is returned. - - Source: Weisstein, Eric W. "Sphere Point Picking." Mathworld. - """ - - theta, u = np.random.uniform(0.0, 2*np.pi, size), \ - np.random.uniform(-1.0, 1.0, size) - - c = np.sqrt(1-u**2) - - if size is None: - return np.array([c*np.cos(theta), c*np.sin(theta), u]) - - points = np.empty((size, 3), dtype) - - points[:,0] = c*np.cos(theta) - points[:,1] = c*np.sin(theta) - points[:,2] = u - - return points diff --git a/photon_map.py b/photon_map.py new file mode 100644 index 0000000..cc07db0 --- /dev/null +++ b/photon_map.py @@ -0,0 +1,86 @@ +import numpy as np +from solid import Solid +from geometry import Geometry +from mesh import mesh_from_stl +from sample import uniform_sphere +from pycuda import autoinit +from pycuda.compiler import SourceModule +from pycuda import gpuarray +import pycuda.driver as cuda +import src +import models +from materials import * + +print 'device %s' % autoinit.device.name() + +lens_mesh = mesh_from_stl(models.dir + '/lens.stl') +lens_solid = Solid(0, lens_mesh, material1=pmt_glass, material2=pmt_vacuum) + +film_mesh = mesh_from_stl(models.dir + '/film.stl') +film_solid = Solid(0, film_mesh, displacement=(-5,0,0), material1=h2o, material2=pmt_vacuum) + +sphere_mesh = mesh_from_stl(models.dir + '/sphere.stl') +sphere_solid = Solid(0, sphere_mesh, displacement=(10,0,0), material1=h2o, material2=pmt_vacuum, surface=aluminum) + +geometry = Geometry() +geometry.add_solid(lens_solid) +geometry.add_solid(film_solid) +geometry.add_solid(sphere_solid) +geometry.build(bits=8) + +module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) +geometry.load(module) +cuda_propagate = module.get_function('propagate') + +nphotons = 100 + +positions = np.tile((10,5,0), nphotons).reshape(nphotons,3) +positions_float3 = np.empty(positions.shape[0], dtype=gpuarray.vec.float3) +positions_float3['x'] = positions[:,0] +positions_float3['y'] = positions[:,1] +positions_float3['z'] = positions[:,2] +positions_gpu = cuda.to_device(positions_float3) + +directions = uniform_sphere(nphotons) +directions_float3 = np.empty(directions.shape[0], dtype=gpuarray.vec.float3) +directions_float3['x'] = directions[:,0] +directions_float3['y'] = directions[:,1] +directions_float3['z'] = directions[:,2] +directions_gpu = cuda.to_device(directions_float3) + +wavelengths = np.random.uniform(200, 800, nphotons).astype(np.float32) +wavelengths_gpu = cuda.to_device(wavelengths) + +times = np.tile(0, nphotons).astype(np.float32) +times_gpu = cuda.to_device(times) + +polarizations = uniform_sphere(nphotons) +polarizations_float3 = np.empty(polarizations.shape[0], dtype=gpuarray.vec.float3) +polarizations_float3['x'] = polarizations[:,0] +polarizations_float3['y'] = polarizations[:,1] +polarizations_float3['z'] = polarizations[:,2] +polarizations_gpu = cuda.to_device(polarizations_float3) + +last_hit_triangles = np.empty(nphotons, dtype=np.int32) +last_hit_triangles_gpu = cuda.to_device(last_hit_triangles) + +states = np.empty(nphotons, dtype=np.int32) +states_gpu = cuda.to_device(states) + +nblocks = 64 + +gpu_kwargs = {'block': (nblocks,1,1), 'grid': (nphotons/nblocks+1,1)} + +cuda_propagate(np.int32(nphotons), positions_gpu, directions_gpu, wavelengths_gpu, polarizations_gpu, times_gpu, states_gpu, last_hit_triangles_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), **gpu_kwargs) + +cuda.memcpy_dtoh(states, states_gpu) + +print states + + + + +#from view import view + +#view(geometry) + diff --git a/sample.py b/sample.py new file mode 100644 index 0000000..b494384 --- /dev/null +++ b/sample.py @@ -0,0 +1,29 @@ +import numpy as np + +def uniform_sphere(size=None, dtype=np.double): + """ + Generate random points isotropically distributed across the unit sphere. + + Args: + - size: int, *optional* + Number of points to generate. If no size is specified, a single + point is returned. + + Source: Weisstein, Eric W. "Sphere Point Picking." Mathworld. + """ + + theta, u = np.random.uniform(0.0, 2*np.pi, size), \ + np.random.uniform(-1.0, 1.0, size) + + c = np.sqrt(1-u**2) + + if size is None: + return np.array([c*np.cos(theta), c*np.sin(theta), u]) + + points = np.empty((size, 3), dtype) + + points[:,0] = c*np.cos(theta) + points[:,1] = c*np.sin(theta) + points[:,2] = u + + return points diff --git a/solid.py b/solid.py index 44f3b3f..929883e 100644 --- a/solid.py +++ b/solid.py @@ -1,9 +1,7 @@ import numpy as np class Solid(object): - def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0), - material1=None, material2=None, surface1=None, surface2=None, - color=0xffffffff): + def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0), material1=None, material2=None, surface=None, color=0xffffffff): self.id = id self.mesh = mesh @@ -33,19 +31,12 @@ class Solid(object): else: self.material2 = np.tile(material2, len(self.mesh)) - if np.iterable(surface1): - if len(surface1) != len(mesh): + if np.iterable(surface): + if len(surface) != len(mesh): raise ValueError('shape mismatch') - self.surface1 = np.array(surface1, dtype=np.object) + self.surface = np.array(surface, dtype=np.object) else: - self.surface1 = np.tile(surface1, len(self.mesh)) - - if np.iterable(surface2): - if len(surface2) != len(mesh): - raise ValueError('shape mismatch') - self.surface2 = np.array(surface2, dtype=np.object) - else: - self.surface2 = np.tile(surface2, len(self.mesh)) + self.surface = np.tile(surface, len(self.mesh)) if np.iterable(color): if len(color) != len(mesh): diff --git a/src/intersect.h b/src/intersect.h index 92bcf0c..e6566f3 100644 --- a/src/intersect.h +++ b/src/intersect.h @@ -1,4 +1,8 @@ //-*-c-*- + +#ifndef __INTERSECT_H__ +#define __INTERSECT_H__ + #include #include "linalg.h" @@ -145,3 +149,5 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con return true; } + +#endif diff --git a/src/kernel.cu b/src/kernel.cu index 796de54..b2e0903 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -7,18 +7,20 @@ #include "rotate.h" #include "intersect.h" #include "materials.h" +#include "physical_constants.h" +#include "random.h" #define STACK_SIZE 500 +#define MAX_DEPTH 5 /* flattened triangle mesh */ texture vertices; __device__ uint4 *triangles; /* material/surface index lookup for each triangle */ -texture material1_index; -texture material2_index; -texture surface1_index; -texture surface2_index; +texture material1_lookup; +texture material2_lookup; +texture surface_lookup; /* lower/upper bounds for the bounding box associated with each node/leaf */ texture upper_bounds; @@ -48,11 +50,10 @@ __device__ bool intersect_node(const float3 &origin, const float3 &direction, co direction `direction` and the global mesh texture. If the ray intersects the texture return the index of the triangle which the ray intersected, else return -1. */ -__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node) +__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node, float &distance) { int triangle_idx = -1; - float distance; float min_distance; if (!intersect_node(origin, direction, start_node)) @@ -133,20 +134,8 @@ __global__ void set_pointer(uint4 *triangle_ptr) /* Initialize random number states */ __global__ void init_rng(unsigned long long seed, unsigned long long offset) { - int idx = blockIdx.x*blockDim.x + threadIdx.x; - curand_init(seed, idx, offset, rng_states+idx); -} - -/* */ -__global__ void uniform_sphere(int nthreads, float3 *points) -{ - int idx = blockIdx.x*blockDim.x + threadIdx.x; - - if (idx >= nthreads) - return; - - //curandState rng = *(rng_states+idx); - + int id = blockIdx.x*blockDim.x + threadIdx.x; + curand_init(seed, id, offset, rng_states+id); } /* Translate `points` by the vector `v` */ @@ -188,7 +177,9 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int float3 direction = *(directions+idx); direction /= norm(direction); - int intersection_idx = intersect_mesh(origin, direction, start_node, first_node); + float distance; + + int intersection_idx = intersect_mesh(origin, direction, start_node, first_node, distance); if (intersection_idx == -1) { @@ -206,23 +197,244 @@ __global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int } } // ray_trace -/* Propagate the photons starting at `origins` traveling in the direction +/* Propagate the photons starting at `positions` traveling in the direction `directions` to their intersection with the global mesh. If the ray intersects the mesh set the hit_solid array value associated with the photon to the triangle index of the triangle the photon intersected, else set the hit_solid array value to -1. */ -__global__ void propagate(int nthreads, float3 *origins, float3 *directions, int start_node, int first_node, int *hit_triangles) +__global__ void propagate(int nthreads, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node) { - int idx = blockIdx.x*blockDim.x + threadIdx.x; + int id = blockIdx.x*blockDim.x + threadIdx.x; - if (idx >= nthreads) + //states[id] = interp_property(materials[0].refractive_index + + //return; + + if (id >= nthreads) return; - float3 origin = *(origins+idx); - float3 direction = *(directions+idx); + float3 position = positions[id]; + float3 direction = directions[id]; direction /= norm(direction); + float3 polarization = polarizations[id]; + float wavelength = wavelengths[id]; + float time = times[id]; + + int last_hit_triangle; + + curandState rng = rng_states[id]; + + unsigned int depth=0; + while (true) + { + depth++; + + if (depth > MAX_DEPTH) + { + states[id] = MAX_DEPTH_REACHED; + break; + } + + float distance; + + last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance); + + if (last_hit_triangle == -1) + { + states[id] = NO_HIT; + break; + } + + uint4 triangle_data = triangles[last_hit_triangle]; + + float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x)); + float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y)); + float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z)); + + int material_in_index = 0xFF & (triangle_data.w >> 24); + int material_out_index = 0xFF & (triangle_data.w >> 16); + int surface_index = 0xFF & (triangle_data.w >> 8); + + if (material_in_index < 0 || material_out_index < 0 || surface_index < 0) + { + states[id] = -2; + break; + } + + float3 surface_normal = cross(v1-v0,v2-v0); + surface_normal /= norm(surface_normal); + + Material material1, material2; + if (dot(surface_normal,-direction) > 0.0f) + { + // outside to inside + material1 = materials[material_out_index]; + material2 = materials[material_in_index]; + } + else + { + // inside to outside + material1 = materials[material_in_index]; + material2 = materials[material_out_index]; + surface_normal = -surface_normal; + } + + float refractive_index1 = interp_property(wavelength, material1.refractive_index); + float refractive_index2 = interp_property(wavelength, material2.refractive_index); + + float absorption_length = interp_property(wavelength, material1.absorption_length); + float scattering_length = interp_property(wavelength, material1.scattering_length); + + float absorption_distance = -absorption_length*logf(curand_uniform(&rng)); + float scattering_distance = -absorption_length*logf(curand_uniform(&rng)); + + if (absorption_distance <= scattering_distance) + { + if (absorption_distance <= distance) + { + time = absorption_distance/(SPEED_OF_LIGHT/refractive_index1); + position = position + absorption_distance*direction; + states[id] = BULK_ABSORB; + break; + } // photon is absorbed in material1 + } + else + { + if (scattering_distance <= distance) + { + time = scattering_distance/(SPEED_OF_LIGHT/refractive_index1); + position = position + scattering_distance*direction; + + float x,y; + while (true) + { + y = curand_uniform(&rng); + x = uniform(0, 2*PI); + + if (y < powf(cosf(x),2)) + break; + } + + float phi = uniform(0, 2*PI); + + float3 old_direction = direction; + float3 old_polarization = polarization; + + direction = rotate(old_direction, x, cross(old_polarization, old_direction)); + polarization = rotate(old_polarization, x, cross(old_polarization, old_direction)); + + direction = rotate(direction, phi, old_polarization); + polarization = rotate(polarization, phi, old_polarization); + + direction /= norm(direction); + polarization /= norm(polarization); + + continue; + } // photon is scattered in material1 + } + + if (surface_index != -1) + { + Surface surface = surfaces[surface_index]; + + float absorption = interp_property(wavelength, surface.absorption); + float reflection_diffuse = interp_property(wavelength, surface.reflection_diffuse); + float reflection_specular = interp_property(wavelength, surface.reflection_specular); + + float sum = absorption + reflection_diffuse + reflection_specular; + absorption /= sum; + reflection_diffuse /= sum; + reflection_specular /= sum; + + float p = curand_uniform(&rng); + + if (p < absorption) + { + // absorb + states[id] = SURFACE_ABSORB; + break; + } + else if (p < absorption + reflection_diffuse) + { + // diffusely reflect + while (true) + { + direction = uniform_sphere(&rng); + + if (dot(direction, surface_normal) > 0.0f) + { + // randomize polarization + polarization = cross(uniform_sphere(&rng), direction); + break; + } + } + + continue; + } + else + { + // specularly reflect + direction = rotate(direction, -PI/2.0f, cross(direction, surface_normal)); + + // polarization ? + + continue; + } + } // surface + + // compute reflection/transmission probability + + // p is normal to the plane of incidence + float3 p = cross(direction, surface_normal); + + float normal_coefficient = dot(polarization, p); + float normal_probability = normal_coefficient*normal_coefficient; + + float incident_angle = acosf(dot(surface_normal,-direction)); + float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2); + + if (curand_uniform(&rng) < normal_probability) + { + // photon polarization normal to plane of incidence + float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle); + float reflection_probability = reflection_coefficient*reflection_coefficient; + + polarization = p; + + if (curand_uniform(&rng) < reflection_probability) + { + // reflect off surface + direction = rotate(surface_normal, PI/2.0f, p); + } + else + { + // transmit + direction = rotate(surface_normal, PI-refracted_angle, p); + } + continue; + } + else + { + // photon polarization parallel to surface + float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle); + float reflection_probability = reflection_coefficient*reflection_coefficient; + + if (curand_uniform(&rng) < reflection_probability) + { + // reflect off surface + direction = rotate(surface_normal, PI/2.0f, p); + polarization = cross(p, direction); + } + else + { + // transmit + direction = rotate(surface_normal, PI-refracted_angle, p); + polarization = cross(p, direction); + } + continue; + } - *(hit_triangles+idx) = intersect_mesh(origin, direction, start_node, first_node); + } // while(true) } // propagate diff --git a/src/materials.h b/src/materials.h index 05aa121..47b7d22 100644 --- a/src/materials.h +++ b/src/materials.h @@ -1,3 +1,20 @@ +#ifndef __MATERIALS_H__ +#define __MATERIALS_H__ + +__device__ float min_wavelength; +__device__ float max_wavelength; +__device__ float wavelength_step; +__device__ unsigned int wavelength_size; + +enum +{ + INIT = -1, + MAX_DEPTH_REACHED, + NO_HIT, + BULK_ABSORB, + SURFACE_ABSORB +}; + struct Material { float *refractive_index; @@ -12,12 +29,34 @@ struct Surface float *reflection_specular; }; -__device__ struct Material materials[100]; -__device__ struct Surface surfaces[100]; +__device__ Material materials[100]; +__device__ Surface surfaces[100]; + +__device__ float interp_property(const float &x, const float *fp) +{ + if (x < min_wavelength) + return fp[0]; + + if (x > max_wavelength) + return fp[wavelength_size-1]; + + unsigned int jl = (x-min_wavelength)/wavelength_step; + float xl = min_wavelength + jl*wavelength_step; + + return xl + (x-xl)*(fp[jl+1]-fp[jl])/wavelength_step; +} extern "C" { +__global__ void set_wavelength_range(float _min_wavelength, float _max_wavelength, float _wavelength_step, unsigned int _wavelength_size) +{ + min_wavelength = _min_wavelength; + max_wavelength = _max_wavelength; + wavelength_step = _wavelength_step; + wavelength_size = _wavelength_size; +} + __global__ void set_material(int material_index, float *refractive_index, float *absorption_length, float *scattering_length) { materials[material_index].refractive_index = refractive_index; @@ -33,3 +72,5 @@ __global__ void set_surface(int surface_index, float *absorption, float *reflect } } // extern "c" + +#endif diff --git a/src/physical_constants.h b/src/physical_constants.h new file mode 100644 index 0000000..2ff87cd --- /dev/null +++ b/src/physical_constants.h @@ -0,0 +1,7 @@ +#ifndef __PHYSICAL_CONSTANTS_H__ +#define __PHYSICAL_CONSTANTS_H__ + +#define SPEED_OF_LIGHT 2.99792458e8 +#define PI 3.141592653589793 + +#endif diff --git a/src/random.h b/src/random.h new file mode 100644 index 0000000..74329da --- /dev/null +++ b/src/random.h @@ -0,0 +1,22 @@ +#ifndef __RANDOM_H__ +#define __RANDOM_H__ + +#include + +#include "physical_constants.h" + +__device__ float uniform(curandState *rng, float a=0.0, float b=1.0) +{ + return a + curand_uniform(rng)*(b-a); +} + +__device__ float3 uniform_sphere(curandState *rng) +{ + float theta = uniform(rng, 0, 2*PI); + float u = uniform(rng, -1, 1); + float c = sqrtf(1-u*u); + + return make_float3(c*cosf(theta), c*sinf(theta), u); +} + +#endif -- cgit