diff --git a/colormaps.c b/colormaps.c new file mode 100644 index 0000000..29ab09e --- /dev/null +++ b/colormaps.c @@ -0,0 +1,1183 @@ +#include "hsluv.c" +#include "turbo_colormap.c" + +/* look-up tables taken from https://github.com/yuki-koyama/tinycolormap */ + +float viridis_srgb_floats[256][3] = { + { 0.267004, 0.004874, 0.329415 }, + { 0.268510, 0.009605, 0.335427 }, + { 0.269944, 0.014625, 0.341379 }, + { 0.271305, 0.019942, 0.347269 }, + { 0.272594, 0.025563, 0.353093 }, + { 0.273809, 0.031497, 0.358853 }, + { 0.274952, 0.037752, 0.364543 }, + { 0.276022, 0.044167, 0.370164 }, + { 0.277018, 0.050344, 0.375715 }, + { 0.277941, 0.056324, 0.381191 }, + { 0.278791, 0.062145, 0.386592 }, + { 0.279566, 0.067836, 0.391917 }, + { 0.280267, 0.073417, 0.397163 }, + { 0.280894, 0.078907, 0.402329 }, + { 0.281446, 0.084320, 0.407414 }, + { 0.281924, 0.089666, 0.412415 }, + { 0.282327, 0.094955, 0.417331 }, + { 0.282656, 0.100196, 0.422160 }, + { 0.282910, 0.105393, 0.426902 }, + { 0.283091, 0.110553, 0.431554 }, + { 0.283197, 0.115680, 0.436115 }, + { 0.283229, 0.120777, 0.440584 }, + { 0.283187, 0.125848, 0.444960 }, + { 0.283072, 0.130895, 0.449241 }, + { 0.282884, 0.135920, 0.453427 }, + { 0.282623, 0.140926, 0.457517 }, + { 0.282290, 0.145912, 0.461510 }, + { 0.281887, 0.150881, 0.465405 }, + { 0.281412, 0.155834, 0.469201 }, + { 0.280868, 0.160771, 0.472899 }, + { 0.280255, 0.165693, 0.476498 }, + { 0.279574, 0.170599, 0.479997 }, + { 0.278826, 0.175490, 0.483397 }, + { 0.278012, 0.180367, 0.486697 }, + { 0.277134, 0.185228, 0.489898 }, + { 0.276194, 0.190074, 0.493001 }, + { 0.275191, 0.194905, 0.496005 }, + { 0.274128, 0.199721, 0.498911 }, + { 0.273006, 0.204520, 0.501721 }, + { 0.271828, 0.209303, 0.504434 }, + { 0.270595, 0.214069, 0.507052 }, + { 0.269308, 0.218818, 0.509577 }, + { 0.267968, 0.223549, 0.512008 }, + { 0.266580, 0.228262, 0.514349 }, + { 0.265145, 0.232956, 0.516599 }, + { 0.263663, 0.237631, 0.518762 }, + { 0.262138, 0.242286, 0.520837 }, + { 0.260571, 0.246922, 0.522828 }, + { 0.258965, 0.251537, 0.524736 }, + { 0.257322, 0.256130, 0.526563 }, + { 0.255645, 0.260703, 0.528312 }, + { 0.253935, 0.265254, 0.529983 }, + { 0.252194, 0.269783, 0.531579 }, + { 0.250425, 0.274290, 0.533103 }, + { 0.248629, 0.278775, 0.534556 }, + { 0.246811, 0.283237, 0.535941 }, + { 0.244972, 0.287675, 0.537260 }, + { 0.243113, 0.292092, 0.538516 }, + { 0.241237, 0.296485, 0.539709 }, + { 0.239346, 0.300855, 0.540844 }, + { 0.237441, 0.305202, 0.541921 }, + { 0.235526, 0.309527, 0.542944 }, + { 0.233603, 0.313828, 0.543914 }, + { 0.231674, 0.318106, 0.544834 }, + { 0.229739, 0.322361, 0.545706 }, + { 0.227802, 0.326594, 0.546532 }, + { 0.225863, 0.330805, 0.547314 }, + { 0.223925, 0.334994, 0.548053 }, + { 0.221989, 0.339161, 0.548752 }, + { 0.220057, 0.343307, 0.549413 }, + { 0.218130, 0.347432, 0.550038 }, + { 0.216210, 0.351535, 0.550627 }, + { 0.214298, 0.355619, 0.551184 }, + { 0.212395, 0.359683, 0.551710 }, + { 0.210503, 0.363727, 0.552206 }, + { 0.208623, 0.367752, 0.552675 }, + { 0.206756, 0.371758, 0.553117 }, + { 0.204903, 0.375746, 0.553533 }, + { 0.203063, 0.379716, 0.553925 }, + { 0.201239, 0.383670, 0.554294 }, + { 0.199430, 0.387607, 0.554642 }, + { 0.197636, 0.391528, 0.554969 }, + { 0.195860, 0.395433, 0.555276 }, + { 0.194100, 0.399323, 0.555565 }, + { 0.192357, 0.403199, 0.555836 }, + { 0.190631, 0.407061, 0.556089 }, + { 0.188923, 0.410910, 0.556326 }, + { 0.187231, 0.414746, 0.556547 }, + { 0.185556, 0.418570, 0.556753 }, + { 0.183898, 0.422383, 0.556944 }, + { 0.182256, 0.426184, 0.557120 }, + { 0.180629, 0.429975, 0.557282 }, + { 0.179019, 0.433756, 0.557430 }, + { 0.177423, 0.437527, 0.557565 }, + { 0.175841, 0.441290, 0.557685 }, + { 0.174274, 0.445044, 0.557792 }, + { 0.172719, 0.448791, 0.557885 }, + { 0.171176, 0.452530, 0.557965 }, + { 0.169646, 0.456262, 0.558030 }, + { 0.168126, 0.459988, 0.558082 }, + { 0.166617, 0.463708, 0.558119 }, + { 0.165117, 0.467423, 0.558141 }, + { 0.163625, 0.471133, 0.558148 }, + { 0.162142, 0.474838, 0.558140 }, + { 0.160665, 0.478540, 0.558115 }, + { 0.159194, 0.482237, 0.558073 }, + { 0.157729, 0.485932, 0.558013 }, + { 0.156270, 0.489624, 0.557936 }, + { 0.154815, 0.493313, 0.557840 }, + { 0.153364, 0.497000, 0.557724 }, + { 0.151918, 0.500685, 0.557587 }, + { 0.150476, 0.504369, 0.557430 }, + { 0.149039, 0.508051, 0.557250 }, + { 0.147607, 0.511733, 0.557049 }, + { 0.146180, 0.515413, 0.556823 }, + { 0.144759, 0.519093, 0.556572 }, + { 0.143343, 0.522773, 0.556295 }, + { 0.141935, 0.526453, 0.555991 }, + { 0.140536, 0.530132, 0.555659 }, + { 0.139147, 0.533812, 0.555298 }, + { 0.137770, 0.537492, 0.554906 }, + { 0.136408, 0.541173, 0.554483 }, + { 0.135066, 0.544853, 0.554029 }, + { 0.133743, 0.548535, 0.553541 }, + { 0.132444, 0.552216, 0.553018 }, + { 0.131172, 0.555899, 0.552459 }, + { 0.129933, 0.559582, 0.551864 }, + { 0.128729, 0.563265, 0.551229 }, + { 0.127568, 0.566949, 0.550556 }, + { 0.126453, 0.570633, 0.549841 }, + { 0.125394, 0.574318, 0.549086 }, + { 0.124395, 0.578002, 0.548287 }, + { 0.123463, 0.581687, 0.547445 }, + { 0.122606, 0.585371, 0.546557 }, + { 0.121831, 0.589055, 0.545623 }, + { 0.121148, 0.592739, 0.544641 }, + { 0.120565, 0.596422, 0.543611 }, + { 0.120092, 0.600104, 0.542530 }, + { 0.119738, 0.603785, 0.541400 }, + { 0.119512, 0.607464, 0.540218 }, + { 0.119423, 0.611141, 0.538982 }, + { 0.119483, 0.614817, 0.537692 }, + { 0.119699, 0.618490, 0.536347 }, + { 0.120081, 0.622161, 0.534946 }, + { 0.120638, 0.625828, 0.533488 }, + { 0.121380, 0.629492, 0.531973 }, + { 0.122312, 0.633153, 0.530398 }, + { 0.123444, 0.636809, 0.528763 }, + { 0.124780, 0.640461, 0.527068 }, + { 0.126326, 0.644107, 0.525311 }, + { 0.128087, 0.647749, 0.523491 }, + { 0.130067, 0.651384, 0.521608 }, + { 0.132268, 0.655014, 0.519661 }, + { 0.134692, 0.658636, 0.517649 }, + { 0.137339, 0.662252, 0.515571 }, + { 0.140210, 0.665859, 0.513427 }, + { 0.143303, 0.669459, 0.511215 }, + { 0.146616, 0.673050, 0.508936 }, + { 0.150148, 0.676631, 0.506589 }, + { 0.153894, 0.680203, 0.504172 }, + { 0.157851, 0.683765, 0.501686 }, + { 0.162016, 0.687316, 0.499129 }, + { 0.166383, 0.690856, 0.496502 }, + { 0.170948, 0.694384, 0.493803 }, + { 0.175707, 0.697900, 0.491033 }, + { 0.180653, 0.701402, 0.488189 }, + { 0.185783, 0.704891, 0.485273 }, + { 0.191090, 0.708366, 0.482284 }, + { 0.196571, 0.711827, 0.479221 }, + { 0.202219, 0.715272, 0.476084 }, + { 0.208030, 0.718701, 0.472873 }, + { 0.214000, 0.722114, 0.469588 }, + { 0.220124, 0.725509, 0.466226 }, + { 0.226397, 0.728888, 0.462789 }, + { 0.232815, 0.732247, 0.459277 }, + { 0.239374, 0.735588, 0.455688 }, + { 0.246070, 0.738910, 0.452024 }, + { 0.252899, 0.742211, 0.448284 }, + { 0.259857, 0.745492, 0.444467 }, + { 0.266941, 0.748751, 0.440573 }, + { 0.274149, 0.751988, 0.436601 }, + { 0.281477, 0.755203, 0.432552 }, + { 0.288921, 0.758394, 0.428426 }, + { 0.296479, 0.761561, 0.424223 }, + { 0.304148, 0.764704, 0.419943 }, + { 0.311925, 0.767822, 0.415586 }, + { 0.319809, 0.770914, 0.411152 }, + { 0.327796, 0.773980, 0.406640 }, + { 0.335885, 0.777018, 0.402049 }, + { 0.344074, 0.780029, 0.397381 }, + { 0.352360, 0.783011, 0.392636 }, + { 0.360741, 0.785964, 0.387814 }, + { 0.369214, 0.788888, 0.382914 }, + { 0.377779, 0.791781, 0.377939 }, + { 0.386433, 0.794644, 0.372886 }, + { 0.395174, 0.797475, 0.367757 }, + { 0.404001, 0.800275, 0.362552 }, + { 0.412913, 0.803041, 0.357269 }, + { 0.421908, 0.805774, 0.351910 }, + { 0.430983, 0.808473, 0.346476 }, + { 0.440137, 0.811138, 0.340967 }, + { 0.449368, 0.813768, 0.335384 }, + { 0.458674, 0.816363, 0.329727 }, + { 0.468053, 0.818921, 0.323998 }, + { 0.477504, 0.821444, 0.318195 }, + { 0.487026, 0.823929, 0.312321 }, + { 0.496615, 0.826376, 0.306377 }, + { 0.506271, 0.828786, 0.300362 }, + { 0.515992, 0.831158, 0.294279 }, + { 0.525776, 0.833491, 0.288127 }, + { 0.535621, 0.835785, 0.281908 }, + { 0.545524, 0.838039, 0.275626 }, + { 0.555484, 0.840254, 0.269281 }, + { 0.565498, 0.842430, 0.262877 }, + { 0.575563, 0.844566, 0.256415 }, + { 0.585678, 0.846661, 0.249897 }, + { 0.595839, 0.848717, 0.243329 }, + { 0.606045, 0.850733, 0.236712 }, + { 0.616293, 0.852709, 0.230052 }, + { 0.626579, 0.854645, 0.223353 }, + { 0.636902, 0.856542, 0.216620 }, + { 0.647257, 0.858400, 0.209861 }, + { 0.657642, 0.860219, 0.203082 }, + { 0.668054, 0.861999, 0.196293 }, + { 0.678489, 0.863742, 0.189503 }, + { 0.688944, 0.865448, 0.182725 }, + { 0.699415, 0.867117, 0.175971 }, + { 0.709898, 0.868751, 0.169257 }, + { 0.720391, 0.870350, 0.162603 }, + { 0.730889, 0.871916, 0.156029 }, + { 0.741388, 0.873449, 0.149561 }, + { 0.751884, 0.874951, 0.143228 }, + { 0.762373, 0.876424, 0.137064 }, + { 0.772852, 0.877868, 0.131109 }, + { 0.783315, 0.879285, 0.125405 }, + { 0.793760, 0.880678, 0.120005 }, + { 0.804182, 0.882046, 0.114965 }, + { 0.814576, 0.883393, 0.110347 }, + { 0.824940, 0.884720, 0.106217 }, + { 0.835270, 0.886029, 0.102646 }, + { 0.845561, 0.887322, 0.099702 }, + { 0.855810, 0.888601, 0.097452 }, + { 0.866013, 0.889868, 0.095953 }, + { 0.876168, 0.891125, 0.095250 }, + { 0.886271, 0.892374, 0.095374 }, + { 0.896320, 0.893616, 0.096335 }, + { 0.906311, 0.894855, 0.098125 }, + { 0.916242, 0.896091, 0.100717 }, + { 0.926106, 0.897330, 0.104071 }, + { 0.935904, 0.898570, 0.108131 }, + { 0.945636, 0.899815, 0.112838 }, + { 0.955300, 0.901065, 0.118128 }, + { 0.964894, 0.902323, 0.123941 }, + { 0.974417, 0.903590, 0.130215 }, + { 0.983868, 0.904867, 0.136897 }, + { 0.993248, 0.906157, 0.143936 } + }; + +float magma_srgb_floats[256][3] = { + { 0.001462, 0.000466, 0.013866 }, + { 0.002258, 0.001295, 0.018331 }, + { 0.003279, 0.002305, 0.023708 }, + { 0.004512, 0.003490, 0.029965 }, + { 0.005950, 0.004843, 0.037130 }, + { 0.007588, 0.006356, 0.044973 }, + { 0.009426, 0.008022, 0.052844 }, + { 0.011465, 0.009828, 0.060750 }, + { 0.013708, 0.011771, 0.068667 }, + { 0.016156, 0.013840, 0.076603 }, + { 0.018815, 0.016026, 0.084584 }, + { 0.021692, 0.018320, 0.092610 }, + { 0.024792, 0.020715, 0.100676 }, + { 0.028123, 0.023201, 0.108787 }, + { 0.031696, 0.025765, 0.116965 }, + { 0.035520, 0.028397, 0.125209 }, + { 0.039608, 0.031090, 0.133515 }, + { 0.043830, 0.033830, 0.141886 }, + { 0.048062, 0.036607, 0.150327 }, + { 0.052320, 0.039407, 0.158841 }, + { 0.056615, 0.042160, 0.167446 }, + { 0.060949, 0.044794, 0.176129 }, + { 0.065330, 0.047318, 0.184892 }, + { 0.069764, 0.049726, 0.193735 }, + { 0.074257, 0.052017, 0.202660 }, + { 0.078815, 0.054184, 0.211667 }, + { 0.083446, 0.056225, 0.220755 }, + { 0.088155, 0.058133, 0.229922 }, + { 0.092949, 0.059904, 0.239164 }, + { 0.097833, 0.061531, 0.248477 }, + { 0.102815, 0.063010, 0.257854 }, + { 0.107899, 0.064335, 0.267289 }, + { 0.113094, 0.065492, 0.276784 }, + { 0.118405, 0.066479, 0.286321 }, + { 0.123833, 0.067295, 0.295879 }, + { 0.129380, 0.067935, 0.305443 }, + { 0.135053, 0.068391, 0.315000 }, + { 0.140858, 0.068654, 0.324538 }, + { 0.146785, 0.068738, 0.334011 }, + { 0.152839, 0.068637, 0.343404 }, + { 0.159018, 0.068354, 0.352688 }, + { 0.165308, 0.067911, 0.361816 }, + { 0.171713, 0.067305, 0.370771 }, + { 0.178212, 0.066576, 0.379497 }, + { 0.184801, 0.065732, 0.387973 }, + { 0.191460, 0.064818, 0.396152 }, + { 0.198177, 0.063862, 0.404009 }, + { 0.204935, 0.062907, 0.411514 }, + { 0.211718, 0.061992, 0.418647 }, + { 0.218512, 0.061158, 0.425392 }, + { 0.225302, 0.060445, 0.431742 }, + { 0.232077, 0.059889, 0.437695 }, + { 0.238826, 0.059517, 0.443256 }, + { 0.245543, 0.059352, 0.448436 }, + { 0.252220, 0.059415, 0.453248 }, + { 0.258857, 0.059706, 0.457710 }, + { 0.265447, 0.060237, 0.461840 }, + { 0.271994, 0.060994, 0.465660 }, + { 0.278493, 0.061978, 0.469190 }, + { 0.284951, 0.063168, 0.472451 }, + { 0.291366, 0.064553, 0.475462 }, + { 0.297740, 0.066117, 0.478243 }, + { 0.304081, 0.067835, 0.480812 }, + { 0.310382, 0.069702, 0.483186 }, + { 0.316654, 0.071690, 0.485380 }, + { 0.322899, 0.073782, 0.487408 }, + { 0.329114, 0.075972, 0.489287 }, + { 0.335308, 0.078236, 0.491024 }, + { 0.341482, 0.080564, 0.492631 }, + { 0.347636, 0.082946, 0.494121 }, + { 0.353773, 0.085373, 0.495501 }, + { 0.359898, 0.087831, 0.496778 }, + { 0.366012, 0.090314, 0.497960 }, + { 0.372116, 0.092816, 0.499053 }, + { 0.378211, 0.095332, 0.500067 }, + { 0.384299, 0.097855, 0.501002 }, + { 0.390384, 0.100379, 0.501864 }, + { 0.396467, 0.102902, 0.502658 }, + { 0.402548, 0.105420, 0.503386 }, + { 0.408629, 0.107930, 0.504052 }, + { 0.414709, 0.110431, 0.504662 }, + { 0.420791, 0.112920, 0.505215 }, + { 0.426877, 0.115395, 0.505714 }, + { 0.432967, 0.117855, 0.506160 }, + { 0.439062, 0.120298, 0.506555 }, + { 0.445163, 0.122724, 0.506901 }, + { 0.451271, 0.125132, 0.507198 }, + { 0.457386, 0.127522, 0.507448 }, + { 0.463508, 0.129893, 0.507652 }, + { 0.469640, 0.132245, 0.507809 }, + { 0.475780, 0.134577, 0.507921 }, + { 0.481929, 0.136891, 0.507989 }, + { 0.488088, 0.139186, 0.508011 }, + { 0.494258, 0.141462, 0.507988 }, + { 0.500438, 0.143719, 0.507920 }, + { 0.506629, 0.145958, 0.507806 }, + { 0.512831, 0.148179, 0.507648 }, + { 0.519045, 0.150383, 0.507443 }, + { 0.525270, 0.152569, 0.507192 }, + { 0.531507, 0.154739, 0.506895 }, + { 0.537755, 0.156894, 0.506551 }, + { 0.544015, 0.159033, 0.506159 }, + { 0.550287, 0.161158, 0.505719 }, + { 0.556571, 0.163269, 0.505230 }, + { 0.562866, 0.165368, 0.504692 }, + { 0.569172, 0.167454, 0.504105 }, + { 0.575490, 0.169530, 0.503466 }, + { 0.581819, 0.171596, 0.502777 }, + { 0.588158, 0.173652, 0.502035 }, + { 0.594508, 0.175701, 0.501241 }, + { 0.600868, 0.177743, 0.500394 }, + { 0.607238, 0.179779, 0.499492 }, + { 0.613617, 0.181811, 0.498536 }, + { 0.620005, 0.183840, 0.497524 }, + { 0.626401, 0.185867, 0.496456 }, + { 0.632805, 0.187893, 0.495332 }, + { 0.639216, 0.189921, 0.494150 }, + { 0.645633, 0.191952, 0.492910 }, + { 0.652056, 0.193986, 0.491611 }, + { 0.658483, 0.196027, 0.490253 }, + { 0.664915, 0.198075, 0.488836 }, + { 0.671349, 0.200133, 0.487358 }, + { 0.677786, 0.202203, 0.485819 }, + { 0.684224, 0.204286, 0.484219 }, + { 0.690661, 0.206384, 0.482558 }, + { 0.697098, 0.208501, 0.480835 }, + { 0.703532, 0.210638, 0.479049 }, + { 0.709962, 0.212797, 0.477201 }, + { 0.716387, 0.214982, 0.475290 }, + { 0.722805, 0.217194, 0.473316 }, + { 0.729216, 0.219437, 0.471279 }, + { 0.735616, 0.221713, 0.469180 }, + { 0.742004, 0.224025, 0.467018 }, + { 0.748378, 0.226377, 0.464794 }, + { 0.754737, 0.228772, 0.462509 }, + { 0.761077, 0.231214, 0.460162 }, + { 0.767398, 0.233705, 0.457755 }, + { 0.773695, 0.236249, 0.455289 }, + { 0.779968, 0.238851, 0.452765 }, + { 0.786212, 0.241514, 0.450184 }, + { 0.792427, 0.244242, 0.447543 }, + { 0.798608, 0.247040, 0.444848 }, + { 0.804752, 0.249911, 0.442102 }, + { 0.810855, 0.252861, 0.439305 }, + { 0.816914, 0.255895, 0.436461 }, + { 0.822926, 0.259016, 0.433573 }, + { 0.828886, 0.262229, 0.430644 }, + { 0.834791, 0.265540, 0.427671 }, + { 0.840636, 0.268953, 0.424666 }, + { 0.846416, 0.272473, 0.421631 }, + { 0.852126, 0.276106, 0.418573 }, + { 0.857763, 0.279857, 0.415496 }, + { 0.863320, 0.283729, 0.412403 }, + { 0.868793, 0.287728, 0.409303 }, + { 0.874176, 0.291859, 0.406205 }, + { 0.879464, 0.296125, 0.403118 }, + { 0.884651, 0.300530, 0.400047 }, + { 0.889731, 0.305079, 0.397002 }, + { 0.894700, 0.309773, 0.393995 }, + { 0.899552, 0.314616, 0.391037 }, + { 0.904281, 0.319610, 0.388137 }, + { 0.908884, 0.324755, 0.385308 }, + { 0.913354, 0.330052, 0.382563 }, + { 0.917689, 0.335500, 0.379915 }, + { 0.921884, 0.341098, 0.377376 }, + { 0.925937, 0.346844, 0.374959 }, + { 0.929845, 0.352734, 0.372677 }, + { 0.933606, 0.358764, 0.370541 }, + { 0.937221, 0.364929, 0.368567 }, + { 0.940687, 0.371224, 0.366762 }, + { 0.944006, 0.377643, 0.365136 }, + { 0.947180, 0.384178, 0.363701 }, + { 0.950210, 0.390820, 0.362468 }, + { 0.953099, 0.397563, 0.361438 }, + { 0.955849, 0.404400, 0.360619 }, + { 0.958464, 0.411324, 0.360014 }, + { 0.960949, 0.418323, 0.359630 }, + { 0.963310, 0.425390, 0.359469 }, + { 0.965549, 0.432519, 0.359529 }, + { 0.967671, 0.439703, 0.359810 }, + { 0.969680, 0.446936, 0.360311 }, + { 0.971582, 0.454210, 0.361030 }, + { 0.973381, 0.461520, 0.361965 }, + { 0.975082, 0.468861, 0.363111 }, + { 0.976690, 0.476226, 0.364466 }, + { 0.978210, 0.483612, 0.366025 }, + { 0.979645, 0.491014, 0.367783 }, + { 0.981000, 0.498428, 0.369734 }, + { 0.982279, 0.505851, 0.371874 }, + { 0.983485, 0.513280, 0.374198 }, + { 0.984622, 0.520713, 0.376698 }, + { 0.985693, 0.528148, 0.379371 }, + { 0.986700, 0.535582, 0.382210 }, + { 0.987646, 0.543015, 0.385210 }, + { 0.988533, 0.550446, 0.388365 }, + { 0.989363, 0.557873, 0.391671 }, + { 0.990138, 0.565296, 0.395122 }, + { 0.990871, 0.572706, 0.398714 }, + { 0.991558, 0.580107, 0.402441 }, + { 0.992196, 0.587502, 0.406299 }, + { 0.992785, 0.594891, 0.410283 }, + { 0.993326, 0.602275, 0.414390 }, + { 0.993834, 0.609644, 0.418613 }, + { 0.994309, 0.616999, 0.422950 }, + { 0.994738, 0.624350, 0.427397 }, + { 0.995122, 0.631696, 0.431951 }, + { 0.995480, 0.639027, 0.436607 }, + { 0.995810, 0.646344, 0.441361 }, + { 0.996096, 0.653659, 0.446213 }, + { 0.996341, 0.660969, 0.451160 }, + { 0.996580, 0.668256, 0.456192 }, + { 0.996775, 0.675541, 0.461314 }, + { 0.996925, 0.682828, 0.466526 }, + { 0.997077, 0.690088, 0.471811 }, + { 0.997186, 0.697349, 0.477182 }, + { 0.997254, 0.704611, 0.482635 }, + { 0.997325, 0.711848, 0.488154 }, + { 0.997351, 0.719089, 0.493755 }, + { 0.997351, 0.726324, 0.499428 }, + { 0.997341, 0.733545, 0.505167 }, + { 0.997285, 0.740772, 0.510983 }, + { 0.997228, 0.747981, 0.516859 }, + { 0.997138, 0.755190, 0.522806 }, + { 0.997019, 0.762398, 0.528821 }, + { 0.996898, 0.769591, 0.534892 }, + { 0.996727, 0.776795, 0.541039 }, + { 0.996571, 0.783977, 0.547233 }, + { 0.996369, 0.791167, 0.553499 }, + { 0.996162, 0.798348, 0.559820 }, + { 0.995932, 0.805527, 0.566202 }, + { 0.995680, 0.812706, 0.572645 }, + { 0.995424, 0.819875, 0.579140 }, + { 0.995131, 0.827052, 0.585701 }, + { 0.994851, 0.834213, 0.592307 }, + { 0.994524, 0.841387, 0.598983 }, + { 0.994222, 0.848540, 0.605696 }, + { 0.993866, 0.855711, 0.612482 }, + { 0.993545, 0.862859, 0.619299 }, + { 0.993170, 0.870024, 0.626189 }, + { 0.992831, 0.877168, 0.633109 }, + { 0.992440, 0.884330, 0.640099 }, + { 0.992089, 0.891470, 0.647116 }, + { 0.991688, 0.898627, 0.654202 }, + { 0.991332, 0.905763, 0.661309 }, + { 0.990930, 0.912915, 0.668481 }, + { 0.990570, 0.920049, 0.675675 }, + { 0.990175, 0.927196, 0.682926 }, + { 0.989815, 0.934329, 0.690198 }, + { 0.989434, 0.941470, 0.697519 }, + { 0.989077, 0.948604, 0.704863 }, + { 0.988717, 0.955742, 0.712242 }, + { 0.988367, 0.962878, 0.719649 }, + { 0.988033, 0.970012, 0.727077 }, + { 0.987691, 0.977154, 0.734536 }, + { 0.987387, 0.984288, 0.742002 }, + { 0.987053, 0.991438, 0.749504 } + }; + +float inferno_srgb_floats[256][3] = { + { 0.001462, 0.000466, 0.013866 }, + { 0.002267, 0.001270, 0.018570 }, + { 0.003299, 0.002249, 0.024239 }, + { 0.004547, 0.003392, 0.030909 }, + { 0.006006, 0.004692, 0.038558 }, + { 0.007676, 0.006136, 0.046836 }, + { 0.009561, 0.007713, 0.055143 }, + { 0.011663, 0.009417, 0.063460 }, + { 0.013995, 0.011225, 0.071862 }, + { 0.016561, 0.013136, 0.080282 }, + { 0.019373, 0.015133, 0.088767 }, + { 0.022447, 0.017199, 0.097327 }, + { 0.025793, 0.019331, 0.105930 }, + { 0.029432, 0.021503, 0.114621 }, + { 0.033385, 0.023702, 0.123397 }, + { 0.037668, 0.025921, 0.132232 }, + { 0.042253, 0.028139, 0.141141 }, + { 0.046915, 0.030324, 0.150164 }, + { 0.051644, 0.032474, 0.159254 }, + { 0.056449, 0.034569, 0.168414 }, + { 0.061340, 0.036590, 0.177642 }, + { 0.066331, 0.038504, 0.186962 }, + { 0.071429, 0.040294, 0.196354 }, + { 0.076637, 0.041905, 0.205799 }, + { 0.081962, 0.043328, 0.215289 }, + { 0.087411, 0.044556, 0.224813 }, + { 0.092990, 0.045583, 0.234358 }, + { 0.098702, 0.046402, 0.243904 }, + { 0.104551, 0.047008, 0.253430 }, + { 0.110536, 0.047399, 0.262912 }, + { 0.116656, 0.047574, 0.272321 }, + { 0.122908, 0.047536, 0.281624 }, + { 0.129285, 0.047293, 0.290788 }, + { 0.135778, 0.046856, 0.299776 }, + { 0.142378, 0.046242, 0.308553 }, + { 0.149073, 0.045468, 0.317085 }, + { 0.155850, 0.044559, 0.325338 }, + { 0.162689, 0.043554, 0.333277 }, + { 0.169575, 0.042489, 0.340874 }, + { 0.176493, 0.041402, 0.348111 }, + { 0.183429, 0.040329, 0.354971 }, + { 0.190367, 0.039309, 0.361447 }, + { 0.197297, 0.038400, 0.367535 }, + { 0.204209, 0.037632, 0.373238 }, + { 0.211095, 0.037030, 0.378563 }, + { 0.217949, 0.036615, 0.383522 }, + { 0.224763, 0.036405, 0.388129 }, + { 0.231538, 0.036405, 0.392400 }, + { 0.238273, 0.036621, 0.396353 }, + { 0.244967, 0.037055, 0.400007 }, + { 0.251620, 0.037705, 0.403378 }, + { 0.258234, 0.038571, 0.406485 }, + { 0.264810, 0.039647, 0.409345 }, + { 0.271347, 0.040922, 0.411976 }, + { 0.277850, 0.042353, 0.414392 }, + { 0.284321, 0.043933, 0.416608 }, + { 0.290763, 0.045644, 0.418637 }, + { 0.297178, 0.047470, 0.420491 }, + { 0.303568, 0.049396, 0.422182 }, + { 0.309935, 0.051407, 0.423721 }, + { 0.316282, 0.053490, 0.425116 }, + { 0.322610, 0.055634, 0.426377 }, + { 0.328921, 0.057827, 0.427511 }, + { 0.335217, 0.060060, 0.428524 }, + { 0.341500, 0.062325, 0.429425 }, + { 0.347771, 0.064616, 0.430217 }, + { 0.354032, 0.066925, 0.430906 }, + { 0.360284, 0.069247, 0.431497 }, + { 0.366529, 0.071579, 0.431994 }, + { 0.372768, 0.073915, 0.432400 }, + { 0.379001, 0.076253, 0.432719 }, + { 0.385228, 0.078591, 0.432955 }, + { 0.391453, 0.080927, 0.433109 }, + { 0.397674, 0.083257, 0.433183 }, + { 0.403894, 0.085580, 0.433179 }, + { 0.410113, 0.087896, 0.433098 }, + { 0.416331, 0.090203, 0.432943 }, + { 0.422549, 0.092501, 0.432714 }, + { 0.428768, 0.094790, 0.432412 }, + { 0.434987, 0.097069, 0.432039 }, + { 0.441207, 0.099338, 0.431594 }, + { 0.447428, 0.101597, 0.431080 }, + { 0.453651, 0.103848, 0.430498 }, + { 0.459875, 0.106089, 0.429846 }, + { 0.466100, 0.108322, 0.429125 }, + { 0.472328, 0.110547, 0.428334 }, + { 0.478558, 0.112764, 0.427475 }, + { 0.484789, 0.114974, 0.426548 }, + { 0.491022, 0.117179, 0.425552 }, + { 0.497257, 0.119379, 0.424488 }, + { 0.503493, 0.121575, 0.423356 }, + { 0.509730, 0.123769, 0.422156 }, + { 0.515967, 0.125960, 0.420887 }, + { 0.522206, 0.128150, 0.419549 }, + { 0.528444, 0.130341, 0.418142 }, + { 0.534683, 0.132534, 0.416667 }, + { 0.540920, 0.134729, 0.415123 }, + { 0.547157, 0.136929, 0.413511 }, + { 0.553392, 0.139134, 0.411829 }, + { 0.559624, 0.141346, 0.410078 }, + { 0.565854, 0.143567, 0.408258 }, + { 0.572081, 0.145797, 0.406369 }, + { 0.578304, 0.148039, 0.404411 }, + { 0.584521, 0.150294, 0.402385 }, + { 0.590734, 0.152563, 0.400290 }, + { 0.596940, 0.154848, 0.398125 }, + { 0.603139, 0.157151, 0.395891 }, + { 0.609330, 0.159474, 0.393589 }, + { 0.615513, 0.161817, 0.391219 }, + { 0.621685, 0.164184, 0.388781 }, + { 0.627847, 0.166575, 0.386276 }, + { 0.633998, 0.168992, 0.383704 }, + { 0.640135, 0.171438, 0.381065 }, + { 0.646260, 0.173914, 0.378359 }, + { 0.652369, 0.176421, 0.375586 }, + { 0.658463, 0.178962, 0.372748 }, + { 0.664540, 0.181539, 0.369846 }, + { 0.670599, 0.184153, 0.366879 }, + { 0.676638, 0.186807, 0.363849 }, + { 0.682656, 0.189501, 0.360757 }, + { 0.688653, 0.192239, 0.357603 }, + { 0.694627, 0.195021, 0.354388 }, + { 0.700576, 0.197851, 0.351113 }, + { 0.706500, 0.200728, 0.347777 }, + { 0.712396, 0.203656, 0.344383 }, + { 0.718264, 0.206636, 0.340931 }, + { 0.724103, 0.209670, 0.337424 }, + { 0.729909, 0.212759, 0.333861 }, + { 0.735683, 0.215906, 0.330245 }, + { 0.741423, 0.219112, 0.326576 }, + { 0.747127, 0.222378, 0.322856 }, + { 0.752794, 0.225706, 0.319085 }, + { 0.758422, 0.229097, 0.315266 }, + { 0.764010, 0.232554, 0.311399 }, + { 0.769556, 0.236077, 0.307485 }, + { 0.775059, 0.239667, 0.303526 }, + { 0.780517, 0.243327, 0.299523 }, + { 0.785929, 0.247056, 0.295477 }, + { 0.791293, 0.250856, 0.291390 }, + { 0.796607, 0.254728, 0.287264 }, + { 0.801871, 0.258674, 0.283099 }, + { 0.807082, 0.262692, 0.278898 }, + { 0.812239, 0.266786, 0.274661 }, + { 0.817341, 0.270954, 0.270390 }, + { 0.822386, 0.275197, 0.266085 }, + { 0.827372, 0.279517, 0.261750 }, + { 0.832299, 0.283913, 0.257383 }, + { 0.837165, 0.288385, 0.252988 }, + { 0.841969, 0.292933, 0.248564 }, + { 0.846709, 0.297559, 0.244113 }, + { 0.851384, 0.302260, 0.239636 }, + { 0.855992, 0.307038, 0.235133 }, + { 0.860533, 0.311892, 0.230606 }, + { 0.865006, 0.316822, 0.226055 }, + { 0.869409, 0.321827, 0.221482 }, + { 0.873741, 0.326906, 0.216886 }, + { 0.878001, 0.332060, 0.212268 }, + { 0.882188, 0.337287, 0.207628 }, + { 0.886302, 0.342586, 0.202968 }, + { 0.890341, 0.347957, 0.198286 }, + { 0.894305, 0.353399, 0.193584 }, + { 0.898192, 0.358911, 0.188860 }, + { 0.902003, 0.364492, 0.184116 }, + { 0.905735, 0.370140, 0.179350 }, + { 0.909390, 0.375856, 0.174563 }, + { 0.912966, 0.381636, 0.169755 }, + { 0.916462, 0.387481, 0.164924 }, + { 0.919879, 0.393389, 0.160070 }, + { 0.923215, 0.399359, 0.155193 }, + { 0.926470, 0.405389, 0.150292 }, + { 0.929644, 0.411479, 0.145367 }, + { 0.932737, 0.417627, 0.140417 }, + { 0.935747, 0.423831, 0.135440 }, + { 0.938675, 0.430091, 0.130438 }, + { 0.941521, 0.436405, 0.125409 }, + { 0.944285, 0.442772, 0.120354 }, + { 0.946965, 0.449191, 0.115272 }, + { 0.949562, 0.455660, 0.110164 }, + { 0.952075, 0.462178, 0.105031 }, + { 0.954506, 0.468744, 0.099874 }, + { 0.956852, 0.475356, 0.094695 }, + { 0.959114, 0.482014, 0.089499 }, + { 0.961293, 0.488716, 0.084289 }, + { 0.963387, 0.495462, 0.079073 }, + { 0.965397, 0.502249, 0.073859 }, + { 0.967322, 0.509078, 0.068659 }, + { 0.969163, 0.515946, 0.063488 }, + { 0.970919, 0.522853, 0.058367 }, + { 0.972590, 0.529798, 0.053324 }, + { 0.974176, 0.536780, 0.048392 }, + { 0.975677, 0.543798, 0.043618 }, + { 0.977092, 0.550850, 0.039050 }, + { 0.978422, 0.557937, 0.034931 }, + { 0.979666, 0.565057, 0.031409 }, + { 0.980824, 0.572209, 0.028508 }, + { 0.981895, 0.579392, 0.026250 }, + { 0.982881, 0.586606, 0.024661 }, + { 0.983779, 0.593849, 0.023770 }, + { 0.984591, 0.601122, 0.023606 }, + { 0.985315, 0.608422, 0.024202 }, + { 0.985952, 0.615750, 0.025592 }, + { 0.986502, 0.623105, 0.027814 }, + { 0.986964, 0.630485, 0.030908 }, + { 0.987337, 0.637890, 0.034916 }, + { 0.987622, 0.645320, 0.039886 }, + { 0.987819, 0.652773, 0.045581 }, + { 0.987926, 0.660250, 0.051750 }, + { 0.987945, 0.667748, 0.058329 }, + { 0.987874, 0.675267, 0.065257 }, + { 0.987714, 0.682807, 0.072489 }, + { 0.987464, 0.690366, 0.079990 }, + { 0.987124, 0.697944, 0.087731 }, + { 0.986694, 0.705540, 0.095694 }, + { 0.986175, 0.713153, 0.103863 }, + { 0.985566, 0.720782, 0.112229 }, + { 0.984865, 0.728427, 0.120785 }, + { 0.984075, 0.736087, 0.129527 }, + { 0.983196, 0.743758, 0.138453 }, + { 0.982228, 0.751442, 0.147565 }, + { 0.981173, 0.759135, 0.156863 }, + { 0.980032, 0.766837, 0.166353 }, + { 0.978806, 0.774545, 0.176037 }, + { 0.977497, 0.782258, 0.185923 }, + { 0.976108, 0.789974, 0.196018 }, + { 0.974638, 0.797692, 0.206332 }, + { 0.973088, 0.805409, 0.216877 }, + { 0.971468, 0.813122, 0.227658 }, + { 0.969783, 0.820825, 0.238686 }, + { 0.968041, 0.828515, 0.249972 }, + { 0.966243, 0.836191, 0.261534 }, + { 0.964394, 0.843848, 0.273391 }, + { 0.962517, 0.851476, 0.285546 }, + { 0.960626, 0.859069, 0.298010 }, + { 0.958720, 0.866624, 0.310820 }, + { 0.956834, 0.874129, 0.323974 }, + { 0.954997, 0.881569, 0.337475 }, + { 0.953215, 0.888942, 0.351369 }, + { 0.951546, 0.896226, 0.365627 }, + { 0.950018, 0.903409, 0.380271 }, + { 0.948683, 0.910473, 0.395289 }, + { 0.947594, 0.917399, 0.410665 }, + { 0.946809, 0.924168, 0.426373 }, + { 0.946392, 0.930761, 0.442367 }, + { 0.946403, 0.937159, 0.458592 }, + { 0.946903, 0.943348, 0.474970 }, + { 0.947937, 0.949318, 0.491426 }, + { 0.949545, 0.955063, 0.507860 }, + { 0.951740, 0.960587, 0.524203 }, + { 0.954529, 0.965896, 0.540361 }, + { 0.957896, 0.971003, 0.556275 }, + { 0.961812, 0.975924, 0.571925 }, + { 0.966249, 0.980678, 0.587206 }, + { 0.971162, 0.985282, 0.602154 }, + { 0.976511, 0.989753, 0.616760 }, + { 0.982257, 0.994109, 0.631017 }, + { 0.988362, 0.998364, 0.644924 } + }; + +float plasma_srgb_floats[256][3] = { + { 0.050383, 0.029803, 0.527975 }, + { 0.063536, 0.028426, 0.533124 }, + { 0.075353, 0.027206, 0.538007 }, + { 0.086222, 0.026125, 0.542658 }, + { 0.096379, 0.025165, 0.547103 }, + { 0.105980, 0.024309, 0.551368 }, + { 0.115124, 0.023556, 0.555468 }, + { 0.123903, 0.022878, 0.559423 }, + { 0.132381, 0.022258, 0.563250 }, + { 0.140603, 0.021687, 0.566959 }, + { 0.148607, 0.021154, 0.570562 }, + { 0.156421, 0.020651, 0.574065 }, + { 0.164070, 0.020171, 0.577478 }, + { 0.171574, 0.019706, 0.580806 }, + { 0.178950, 0.019252, 0.584054 }, + { 0.186213, 0.018803, 0.587228 }, + { 0.193374, 0.018354, 0.590330 }, + { 0.200445, 0.017902, 0.593364 }, + { 0.207435, 0.017442, 0.596333 }, + { 0.214350, 0.016973, 0.599239 }, + { 0.221197, 0.016497, 0.602083 }, + { 0.227983, 0.016007, 0.604867 }, + { 0.234715, 0.015502, 0.607592 }, + { 0.241396, 0.014979, 0.610259 }, + { 0.248032, 0.014439, 0.612868 }, + { 0.254627, 0.013882, 0.615419 }, + { 0.261183, 0.013308, 0.617911 }, + { 0.267703, 0.012716, 0.620346 }, + { 0.274191, 0.012109, 0.622722 }, + { 0.280648, 0.011488, 0.625038 }, + { 0.287076, 0.010855, 0.627295 }, + { 0.293478, 0.010213, 0.629490 }, + { 0.299855, 0.009561, 0.631624 }, + { 0.306210, 0.008902, 0.633694 }, + { 0.312543, 0.008239, 0.635700 }, + { 0.318856, 0.007576, 0.637640 }, + { 0.325150, 0.006915, 0.639512 }, + { 0.331426, 0.006261, 0.641316 }, + { 0.337683, 0.005618, 0.643049 }, + { 0.343925, 0.004991, 0.644710 }, + { 0.350150, 0.004382, 0.646298 }, + { 0.356359, 0.003798, 0.647810 }, + { 0.362553, 0.003243, 0.649245 }, + { 0.368733, 0.002724, 0.650601 }, + { 0.374897, 0.002245, 0.651876 }, + { 0.381047, 0.001814, 0.653068 }, + { 0.387183, 0.001434, 0.654177 }, + { 0.393304, 0.001114, 0.655199 }, + { 0.399411, 0.000859, 0.656133 }, + { 0.405503, 0.000678, 0.656977 }, + { 0.411580, 0.000577, 0.657730 }, + { 0.417642, 0.000564, 0.658390 }, + { 0.423689, 0.000646, 0.658956 }, + { 0.429719, 0.000831, 0.659425 }, + { 0.435734, 0.001127, 0.659797 }, + { 0.441732, 0.001540, 0.660069 }, + { 0.447714, 0.002080, 0.660240 }, + { 0.453677, 0.002755, 0.660310 }, + { 0.459623, 0.003574, 0.660277 }, + { 0.465550, 0.004545, 0.660139 }, + { 0.471457, 0.005678, 0.659897 }, + { 0.477344, 0.006980, 0.659549 }, + { 0.483210, 0.008460, 0.659095 }, + { 0.489055, 0.010127, 0.658534 }, + { 0.494877, 0.011990, 0.657865 }, + { 0.500678, 0.014055, 0.657088 }, + { 0.506454, 0.016333, 0.656202 }, + { 0.512206, 0.018833, 0.655209 }, + { 0.517933, 0.021563, 0.654109 }, + { 0.523633, 0.024532, 0.652901 }, + { 0.529306, 0.027747, 0.651586 }, + { 0.534952, 0.031217, 0.650165 }, + { 0.540570, 0.034950, 0.648640 }, + { 0.546157, 0.038954, 0.647010 }, + { 0.551715, 0.043136, 0.645277 }, + { 0.557243, 0.047331, 0.643443 }, + { 0.562738, 0.051545, 0.641509 }, + { 0.568201, 0.055778, 0.639477 }, + { 0.573632, 0.060028, 0.637349 }, + { 0.579029, 0.064296, 0.635126 }, + { 0.584391, 0.068579, 0.632812 }, + { 0.589719, 0.072878, 0.630408 }, + { 0.595011, 0.077190, 0.627917 }, + { 0.600266, 0.081516, 0.625342 }, + { 0.605485, 0.085854, 0.622686 }, + { 0.610667, 0.090204, 0.619951 }, + { 0.615812, 0.094564, 0.617140 }, + { 0.620919, 0.098934, 0.614257 }, + { 0.625987, 0.103312, 0.611305 }, + { 0.631017, 0.107699, 0.608287 }, + { 0.636008, 0.112092, 0.605205 }, + { 0.640959, 0.116492, 0.602065 }, + { 0.645872, 0.120898, 0.598867 }, + { 0.650746, 0.125309, 0.595617 }, + { 0.655580, 0.129725, 0.592317 }, + { 0.660374, 0.134144, 0.588971 }, + { 0.665129, 0.138566, 0.585582 }, + { 0.669845, 0.142992, 0.582154 }, + { 0.674522, 0.147419, 0.578688 }, + { 0.679160, 0.151848, 0.575189 }, + { 0.683758, 0.156278, 0.571660 }, + { 0.688318, 0.160709, 0.568103 }, + { 0.692840, 0.165141, 0.564522 }, + { 0.697324, 0.169573, 0.560919 }, + { 0.701769, 0.174005, 0.557296 }, + { 0.706178, 0.178437, 0.553657 }, + { 0.710549, 0.182868, 0.550004 }, + { 0.714883, 0.187299, 0.546338 }, + { 0.719181, 0.191729, 0.542663 }, + { 0.723444, 0.196158, 0.538981 }, + { 0.727670, 0.200586, 0.535293 }, + { 0.731862, 0.205013, 0.531601 }, + { 0.736019, 0.209439, 0.527908 }, + { 0.740143, 0.213864, 0.524216 }, + { 0.744232, 0.218288, 0.520524 }, + { 0.748289, 0.222711, 0.516834 }, + { 0.752312, 0.227133, 0.513149 }, + { 0.756304, 0.231555, 0.509468 }, + { 0.760264, 0.235976, 0.505794 }, + { 0.764193, 0.240396, 0.502126 }, + { 0.768090, 0.244817, 0.498465 }, + { 0.771958, 0.249237, 0.494813 }, + { 0.775796, 0.253658, 0.491171 }, + { 0.779604, 0.258078, 0.487539 }, + { 0.783383, 0.262500, 0.483918 }, + { 0.787133, 0.266922, 0.480307 }, + { 0.790855, 0.271345, 0.476706 }, + { 0.794549, 0.275770, 0.473117 }, + { 0.798216, 0.280197, 0.469538 }, + { 0.801855, 0.284626, 0.465971 }, + { 0.805467, 0.289057, 0.462415 }, + { 0.809052, 0.293491, 0.458870 }, + { 0.812612, 0.297928, 0.455338 }, + { 0.816144, 0.302368, 0.451816 }, + { 0.819651, 0.306812, 0.448306 }, + { 0.823132, 0.311261, 0.444806 }, + { 0.826588, 0.315714, 0.441316 }, + { 0.830018, 0.320172, 0.437836 }, + { 0.833422, 0.324635, 0.434366 }, + { 0.836801, 0.329105, 0.430905 }, + { 0.840155, 0.333580, 0.427455 }, + { 0.843484, 0.338062, 0.424013 }, + { 0.846788, 0.342551, 0.420579 }, + { 0.850066, 0.347048, 0.417153 }, + { 0.853319, 0.351553, 0.413734 }, + { 0.856547, 0.356066, 0.410322 }, + { 0.859750, 0.360588, 0.406917 }, + { 0.862927, 0.365119, 0.403519 }, + { 0.866078, 0.369660, 0.400126 }, + { 0.869203, 0.374212, 0.396738 }, + { 0.872303, 0.378774, 0.393355 }, + { 0.875376, 0.383347, 0.389976 }, + { 0.878423, 0.387932, 0.386600 }, + { 0.881443, 0.392529, 0.383229 }, + { 0.884436, 0.397139, 0.379860 }, + { 0.887402, 0.401762, 0.376494 }, + { 0.890340, 0.406398, 0.373130 }, + { 0.893250, 0.411048, 0.369768 }, + { 0.896131, 0.415712, 0.366407 }, + { 0.898984, 0.420392, 0.363047 }, + { 0.901807, 0.425087, 0.359688 }, + { 0.904601, 0.429797, 0.356329 }, + { 0.907365, 0.434524, 0.352970 }, + { 0.910098, 0.439268, 0.349610 }, + { 0.912800, 0.444029, 0.346251 }, + { 0.915471, 0.448807, 0.342890 }, + { 0.918109, 0.453603, 0.339529 }, + { 0.920714, 0.458417, 0.336166 }, + { 0.923287, 0.463251, 0.332801 }, + { 0.925825, 0.468103, 0.329435 }, + { 0.928329, 0.472975, 0.326067 }, + { 0.930798, 0.477867, 0.322697 }, + { 0.933232, 0.482780, 0.319325 }, + { 0.935630, 0.487712, 0.315952 }, + { 0.937990, 0.492667, 0.312575 }, + { 0.940313, 0.497642, 0.309197 }, + { 0.942598, 0.502639, 0.305816 }, + { 0.944844, 0.507658, 0.302433 }, + { 0.947051, 0.512699, 0.299049 }, + { 0.949217, 0.517763, 0.295662 }, + { 0.951344, 0.522850, 0.292275 }, + { 0.953428, 0.527960, 0.288883 }, + { 0.955470, 0.533093, 0.285490 }, + { 0.957469, 0.538250, 0.282096 }, + { 0.959424, 0.543431, 0.278701 }, + { 0.961336, 0.548636, 0.275305 }, + { 0.963203, 0.553865, 0.271909 }, + { 0.965024, 0.559118, 0.268513 }, + { 0.966798, 0.564396, 0.265118 }, + { 0.968526, 0.569700, 0.261721 }, + { 0.970205, 0.575028, 0.258325 }, + { 0.971835, 0.580382, 0.254931 }, + { 0.973416, 0.585761, 0.251540 }, + { 0.974947, 0.591165, 0.248151 }, + { 0.976428, 0.596595, 0.244767 }, + { 0.977856, 0.602051, 0.241387 }, + { 0.979233, 0.607532, 0.238013 }, + { 0.980556, 0.613039, 0.234646 }, + { 0.981826, 0.618572, 0.231287 }, + { 0.983041, 0.624131, 0.227937 }, + { 0.984199, 0.629718, 0.224595 }, + { 0.985301, 0.635330, 0.221265 }, + { 0.986345, 0.640969, 0.217948 }, + { 0.987332, 0.646633, 0.214648 }, + { 0.988260, 0.652325, 0.211364 }, + { 0.989128, 0.658043, 0.208100 }, + { 0.989935, 0.663787, 0.204859 }, + { 0.990681, 0.669558, 0.201642 }, + { 0.991365, 0.675355, 0.198453 }, + { 0.991985, 0.681179, 0.195295 }, + { 0.992541, 0.687030, 0.192170 }, + { 0.993032, 0.692907, 0.189084 }, + { 0.993456, 0.698810, 0.186041 }, + { 0.993814, 0.704741, 0.183043 }, + { 0.994103, 0.710698, 0.180097 }, + { 0.994324, 0.716681, 0.177208 }, + { 0.994474, 0.722691, 0.174381 }, + { 0.994553, 0.728728, 0.171622 }, + { 0.994561, 0.734791, 0.168938 }, + { 0.994495, 0.740880, 0.166335 }, + { 0.994355, 0.746995, 0.163821 }, + { 0.994141, 0.753137, 0.161404 }, + { 0.993851, 0.759304, 0.159092 }, + { 0.993482, 0.765499, 0.156891 }, + { 0.993033, 0.771720, 0.154808 }, + { 0.992505, 0.777967, 0.152855 }, + { 0.991897, 0.784239, 0.151042 }, + { 0.991209, 0.790537, 0.149377 }, + { 0.990439, 0.796859, 0.147870 }, + { 0.989587, 0.803205, 0.146529 }, + { 0.988648, 0.809579, 0.145357 }, + { 0.987621, 0.815978, 0.144363 }, + { 0.986509, 0.822401, 0.143557 }, + { 0.985314, 0.828846, 0.142945 }, + { 0.984031, 0.835315, 0.142528 }, + { 0.982653, 0.841812, 0.142303 }, + { 0.981190, 0.848329, 0.142279 }, + { 0.979644, 0.854866, 0.142453 }, + { 0.977995, 0.861432, 0.142808 }, + { 0.976265, 0.868016, 0.143351 }, + { 0.974443, 0.874622, 0.144061 }, + { 0.972530, 0.881250, 0.144923 }, + { 0.970533, 0.887896, 0.145919 }, + { 0.968443, 0.894564, 0.147014 }, + { 0.966271, 0.901249, 0.148180 }, + { 0.964021, 0.907950, 0.149370 }, + { 0.961681, 0.914672, 0.150520 }, + { 0.959276, 0.921407, 0.151566 }, + { 0.956808, 0.928152, 0.152409 }, + { 0.954287, 0.934908, 0.152921 }, + { 0.951726, 0.941671, 0.152925 }, + { 0.949151, 0.948435, 0.152178 }, + { 0.946602, 0.955190, 0.150328 }, + { 0.944152, 0.961916, 0.146861 }, + { 0.941896, 0.968590, 0.140956 }, + { 0.940015, 0.975158, 0.131326 } + }; + + +void hsl_to_rgb_jet(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */ +/* h = hue, s = saturation, l = luminosity */ +{ + double c = 0.0, m = 0.0, x = 0.0; + + c = (1.0 - fabs(2.0 * l - 1.0)) * s; + m = 1.0 * (l - 0.5 * c); + x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0)); + + if (h >= 0.0 && h < 60.0) + { + rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m; + } + else if (h < 120.0) + { + rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m; + } + else if (h < 180.0) + { + rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m; + } + else if (h < 240.0) + { + rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m; + } + else if (h < 300.0) + { + rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m; + } + else if (h < 360.0) + { + rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m; + } + else + { + rgb[0] = m; rgb[1] = m; rgb[2] = m; + } +} + +void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */ +{ + int color_hue, i; + double r, g, b, ratio = 0.711111111; /* ratio equals 256/360 */ + + switch (COLOR_PALETTE) { + case (COL_JET): + { + hsl_to_rgb_jet(h, s, l, rgb); + break; + } + case (COL_HSLUV): + { + hsluv2rgb(h, 100.0*s, l, &r, &g, &b); + rgb[0] = 2.0*l*r*100.0; + rgb[1] = 2.0*l*g*100.0; + rgb[2] = 2.0*l*b*100.0; + break; + } + case (COL_TURBO): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)turbo_srgb_floats[color_hue][i]; + break; + } + case (COL_VIRIDIS): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)viridis_srgb_floats[color_hue][i]; + break; + } + case (COL_MAGMA): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)magma_srgb_floats[color_hue][i]; + break; + } + case (COL_INFERNO): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)inferno_srgb_floats[color_hue][i]; + break; + } + case (COL_PLASMA): + { + color_hue = 255 - (int)(ratio*h); + for (i=0; i<3; i++) rgb[i] = 2.0*l*(double)plasma_srgb_floats[color_hue][i]; + break; + } + } +} + +void amp_to_rgb(double h, double rgb[3]) /* color conversion from amplitude in [0,1] to RGB */ +{ + int color_hue, i; + double r, g, b; + + color_hue = (int)(256.0*h); + if (color_hue > 255) color_hue = 255; + + switch (COLOR_PALETTE) { + case (COL_JET): + { + hsl_to_rgb_jet(360.0*(1.0 - h), 0.9, 0.5, rgb); + break; + } + case (COL_HSLUV): + { + hsluv2rgb(360.0*(1.0 - h), 90.0, 0.5, &r, &g, &b); + rgb[0] = r*100.0; + rgb[1] = g*100.0; + rgb[2] = b*100.0; + break; + } + case (COL_TURBO): + { + for (i=0; i<3; i++) rgb[i] = (double)turbo_srgb_floats[color_hue][i]; + break; + } + case (COL_VIRIDIS): + { + for (i=0; i<3; i++) rgb[i] = (double)viridis_srgb_floats[color_hue][i]; + break; + } + case (COL_MAGMA): + { + for (i=0; i<3; i++) rgb[i] = (double)magma_srgb_floats[color_hue][i]; + break; + } + case (COL_INFERNO): + { + for (i=0; i<3; i++) rgb[i] = (double)inferno_srgb_floats[color_hue][i]; + break; + } + case (COL_PLASMA): + { + for (i=0; i<3; i++) rgb[i] = (double)plasma_srgb_floats[color_hue][i]; + break; + } + } +} + diff --git a/colors_waves.c b/colors_waves.c new file mode 100644 index 0000000..af101b2 --- /dev/null +++ b/colors_waves.c @@ -0,0 +1,130 @@ +#include "colormaps.c" + +double color_amplitude(double value, double scale, int time) +/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */ +{ + return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION))); +} + +double color_amplitude_asym(double value, double scale, int time) +/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */ +{ + return(2.0*tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)) - 1.0); +} + + +void color_scheme(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */ +{ + double hue, y, r, amplitude; + int intpart; + + /* saturation = r, luminosity = y */ + switch (scheme) { + case (C_LUM): + { + hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = LUMMEAN + amplitude*LUMAMP; + intpart = (int)y; + y -= (double)intpart; + hsl_to_rgb(hue, r, y, rgb); + break; + } + case (C_HUE): + { + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = 0.5; + hue = HUEMEAN + amplitude*HUEAMP; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + hsl_to_rgb(hue, r, y, rgb); + break; + } + case (C_ONEDIM): + { + amplitude = color_amplitude(value, scale, time); + amp_to_rgb(0.5*(1.0 + amplitude), rgb); + break; + } + } +} + + +void color_scheme_lum(int scheme, double value, double scale, int time, double lum, double rgb[3]) /* color scheme */ +{ + double hue, y, r, amplitude; + int intpart; + + /* saturation = r, luminosity = y */ + switch (scheme) { + case (C_LUM): + { + hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = LUMMEAN + amplitude*LUMAMP; + intpart = (int)y; + y -= (double)intpart; + hsl_to_rgb(hue, r, y, rgb); + break; + } + case (C_HUE): + { + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = lum; + hue = HUEMEAN + amplitude*HUEAMP; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + hsl_to_rgb(hue, r, y, rgb); + break; + } + } +} + +void color_scheme_asym(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */ +{ + double hue, y, r, amplitude; + int intpart; + + /* saturation = r, luminosity = y */ + switch (scheme) { + case (C_LUM): + { + hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + r = 0.9; + amplitude = color_amplitude(value, scale, time); + y = LUMMEAN + amplitude*LUMAMP; + intpart = (int)y; + y -= (double)intpart; + hsl_to_rgb(hue, r, y, rgb); + break; + } + case (C_HUE): + { + r = 0.9; + amplitude = color_amplitude_asym(value, scale, time); + y = 0.5; + hue = HUEMEAN + 0.8*amplitude*HUEAMP; + if (hue < 0.0) hue += 360.0; + if (hue >= 360.0) hue -= 360.0; + hsl_to_rgb(hue, r, y, rgb); + break; + } + case (C_ONEDIM): + { + amplitude = color_amplitude(value, scale, time); + amp_to_rgb(amplitude, rgb); + break; + } + } +} + diff --git a/drop_billiard.c b/drop_billiard.c index ef8741c..883f7a1 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -30,14 +30,14 @@ #define MOVIE 0 /* set to 1 to generate movie */ -// #define WINWIDTH 720 /* window width */ -#define WINWIDTH 1280 /* window width */ +#define WINWIDTH 720 /* window width */ +// #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -// #define XMIN -1.125 -// #define XMAX 1.125 /* x interval */ +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +#define XMIN -1.125 +#define XMAX 1.125 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ @@ -65,9 +65,9 @@ // #define LAMBDA 3.75738973 /* sin(36°)/sin(9°) for 5-star shape with 90° angles */ // #define LAMBDA -1.73205080756888 /* -sqrt(3) for Reuleaux triangle */ // #define LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */ -#define MU 0.9 /* second parameter controlling shape of billiard */ +#define MU 1.0 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NPOLY 4 /* number of sides of polygon */ +#define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ #define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */ @@ -78,10 +78,10 @@ #define NPART 10000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ -#define NSTEPS 4200 /* number of frames of movie */ -#define TIME 100 /* time between movie frames, for fluidity of real-time simulation */ +#define NSTEPS 5000 /* number of frames of movie */ +#define TIME 150 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.0001 /* integration step */ -#define NVID 50 /* number of iterations between images displayed on screen */ +#define NVID 75 /* number of iterations between images displayed on screen */ /* Decreasing TIME accelerates the animation and the movie */ /* For constant speed of movie, TIME*DPHI should be kept constant */ @@ -99,12 +99,14 @@ /* color and other graphical parameters */ -#define NCOLORS 12 /* number of colors */ -#define COLORSHIFT 0 /* hue of initial color */ +#define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */ + +#define NCOLORS 6 /* number of colors */ +#define COLORSHIFT 3 /* hue of initial color */ #define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define NSEG 100 /* number of segments of boundary */ #define BILLIARD_WIDTH 4 /* width of billiard */ -#define FRONT_WIDTH 3 /* width of wave front */ +#define FRONT_WIDTH 4 /* width of wave front */ #define BLACK 1 /* set to 1 for black background */ #define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ diff --git a/global_particles.c b/global_particles.c index 7ac20f2..8e4cbfa 100644 --- a/global_particles.c +++ b/global_particles.c @@ -36,6 +36,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define D_CIRCLES 20 /* several circles */ #define D_CIRCLES_IN_RECT 21 /* several circles inside a rectangle */ #define D_CIRCLES_IN_GENUSN 22 /* several circles in polygon with identified opposite sides */ +#define D_CIRCLES_IN_TORUS 23 /* several circles in a rectangle with periodic boundary conditions */ #define C_FOUR_CIRCLES 0 /* four circles almost touching each other */ #define C_SQUARE 1 /* square grid of circles */ @@ -49,4 +50,17 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define C_LASER 11 /* laser fight in a room of mirrors */ -#define C_LASER_GENUSN 12 /* laser fight in a translation surface */ \ No newline at end of file +#define C_LASER_GENUSN 12 /* laser fight in a translation surface */ + +/* Color palettes */ + +#define COL_JET 0 /* JET color palette */ +#define COL_HSLUV 1 /* HSLUV color palette (perceptually uniform) */ + +#define COL_TURBO 10 /* TURBO color palette (by Anton Mikhailov) */ +#define COL_VIRIDIS 11 /* Viridis color palette */ +#define COL_MAGMA 12 /* Magma color palette */ +#define COL_INFERNO 13 /* Inferno color palette */ +#define COL_PLASMA 14 /* Plasma color palette */ + + diff --git a/global_pdes.c b/global_pdes.c index 5baca82..db627eb 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -1,5 +1,8 @@ /* Global variables and parameters for wave_billiard, heat and schrodinger */ +// #include "hsluv.c" + + /* Basic math */ #define PI 3.141592654 @@ -33,13 +36,17 @@ #define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */ #define D_CIRCLES 20 /* several circles */ +#define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ #define D_FOUR_PARABOLAS 31 /* four parabolas with axes in NSEW directions */ #define D_POLY_PARABOLAS 32 /* polygon with parabolic sides */ -#define D_PENROSE 33 /* Penrose illumination problem */ +#define D_PENROSE 33 /* Penrose unilluminable room */ #define D_HYPERBOLA 34 /* one branch of hyperbola */ +#define D_TOKARSKY 35 /* Tokarsky unilluminable room */ -#define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */ +#define D_ISOSPECTRAL 37 /* isospectral billiards */ + +#define NMAXCIRCLES 600 /* total number of circles (must be at least NCX*NCY for square grid) */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -48,7 +55,7 @@ #define C_RAND_POISSON 4 /* random Poisson point process */ #define C_CLOAK 5 /* invisibility cloak */ #define C_CLOAK_A 6 /* first optimized invisibility cloak */ - +#define C_LASER 7 /* laser fight in a room of mirrors */ #define C_POISSON_DISC 8 /* Poisson disc sampling */ #define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */ @@ -102,6 +109,18 @@ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ #define C_HUE 1 /* color scheme modifies hue */ #define C_PHASE 2 /* color scheme shows phase */ +#define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */ + +/* Color palettes */ + +#define COL_JET 0 /* JET color palette */ +#define COL_HSLUV 1 /* HSLUV color palette (perceptually uniform) */ + +#define COL_TURBO 10 /* TURBO color palette (by Anton Mikhailov) */ +#define COL_VIRIDIS 11 /* Viridis color palette */ +#define COL_MAGMA 12 /* Magma color palette */ +#define COL_INFERNO 13 /* Inferno color palette */ +#define COL_PLASMA 14 /* Plasma color palette */ double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES]; /* position and radius of circular scatterers */ diff --git a/heat.c b/heat.c index 9cf37a6..b05b4bd 100644 --- a/heat.c +++ b/heat.c @@ -80,6 +80,17 @@ #define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard in sub_wave.c */ @@ -153,11 +164,10 @@ // #define HUEAMP -130.0 /* amplitude of variation of hue for color scheme C_HUE */ -#include "hsluv.c" - #include "global_pdes.c" #include "sub_wave.c" + double courant2; /* Courant parameter squared */ double dx2; /* spatial step size squared */ double intstep; /* integration step */ diff --git a/mangrove.c b/mangrove.c index f3bc8d1..74c2632 100644 --- a/mangrove.c +++ b/mangrove.c @@ -49,8 +49,10 @@ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -#define NX 1280 /* number of grid points on x axis */ -#define NY 720 /* number of grid points on y axis */ +#define NX 640 /* number of grid points on x axis */ +#define NY 360 /* number of grid points on y axis */ +// #define NX 1280 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -63,6 +65,7 @@ #define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +// #define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ @@ -80,6 +83,17 @@ #define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ @@ -90,18 +104,24 @@ #define OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define X_SHIFT -0.9 /* x range on which to apply OSCILLATE_TOPBOT */ -#define OMEGA 0.002 /* frequency of periodic excitation */ -#define K_BC 3.0 /* spatial period of periodic excitation in y direction */ -#define KX_BC 30.0 /* spatial period of periodic excitation in x direction */ -#define KY_BC 10.0 /* spatial period of periodic excitation in y direction */ +#define OMEGA 0.00133333333 /* frequency of periodic excitation */ +#define K_BC 3.0 /* spatial period of periodic excitation in y direction */ +#define KX_BC 10.0 /* spatial period of periodic excitation in x direction */ +#define KY_BC 3.3333 /* spatial period of periodic excitation in y direction */ +// #define KX_BC 20.0 /* spatial period of periodic excitation in x direction */ +// #define KY_BC 6.66666 /* spatial period of periodic excitation in y direction */ +// #define OMEGA 0.002 /* frequency of periodic excitation */ +// #define K_BC 3.0 /* spatial period of periodic excitation in y direction */ +// #define KX_BC 30.0 /* spatial period of periodic excitation in x direction */ +// #define KY_BC 10.0 /* spatial period of periodic excitation in y direction */ #define AMPLITUDE 1.0 /* amplitude of periodic excitation */ #define COURANT 0.02 /* Courant number */ -#define COURANTB 0.01 /* Courant number in medium B */ +#define COURANTB 0.015 /* Courant number in medium B */ // #define COURANTB 0.00666 /* Courant number in medium B */ -#define GAMMA 2.0e-6 /* damping factor in wave equation */ -#define GAMMAB 2.5e-4 /* damping factor in wave equation */ -// #define GAMMAB 5.0e-4 /* damping factor in wave equation */ -// #define GAMMAB 1.0e-4 /* damping factor in wave equation */ +#define GAMMA 3.0e-6 /* damping factor in wave equation */ +#define GAMMAB 5.0e-4 /* damping factor in wave equation */ +// #define GAMMA 2.0e-6 /* damping factor in wave equation */ +// #define GAMMAB 2.5e-4 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -119,8 +139,8 @@ /* Parameters for length and speed of simulation */ -// #define NSTEPS 1000 /* number of frames of movie */ -#define NSTEPS 4500 /* number of frames of movie */ +#define NSTEPS 1000 /* number of frames of movie */ +// #define NSTEPS 5500 /* number of frames of movie */ #define NVID 60 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 100 /* time after which to start saving frames */ @@ -167,42 +187,40 @@ #define MANGROVE_HUE_MIN 180.0 /* color of original mangrove */ #define MANGROVE_HUE_MAX -50.0 /* color of saturated mangrove */ // #define MANGROVE_EMAX 5.0e-3 /* max energy for mangrove to survive */ -#define MANGROVE_EMAX 1.1e-3 /* max energy for mangrove to survive */ +#define MANGROVE_EMAX 1.5e-3 /* max energy for mangrove to survive */ #define RANDOM_RADIUS 1 /* set to 1 for random circle radius */ -#define ERODE_MANGROVES 0 /* set to 1 for mangroves to be eroded */ +#define ERODE_MANGROVES 1 /* set to 1 for mangroves to be eroded */ +#define RECOVER_MANGROVES 1 /* set to 1 to allow mangroves to recover */ #define MOVE_MANGROVES 1 /* set to 1 for mobile mangroves */ #define DETACH_MANGROVES 1 /* set to 1 for mangroves to be able to detach */ #define INERTIA 1 /* set to 1 for taking inertia into account */ +#define REPELL_MANGROVES 1 /* set to 1 for mangroves to repell each other */ #define DT_MANGROVE 0.1 /* time step for mangrove displacement */ -#define KSPRING 0.25 /* spring constant of mangroves */ -#define KWAVE 2.0 /* constant in force due to wave gradient */ +#define KSPRING 0.05 /* spring constant of mangroves */ +#define KWAVE 4.0 /* constant in force due to wave gradient */ +#define KREPEL 5.0 /* constant in repelling force between mangroves */ +#define REPEL_RADIUS 1.1 /* radius in which repelling force acts (in units of mangrove radius) */ #define DXMAX 0.02 /* max displacement of mangrove in one time step */ -#define L_DETACH 0.2 /* spring length beyond which mangroves detach */ -#define DAMP_MANGROVE 0.2 /* damping coefficient of mangroves */ +#define L_DETACH 0.25 /* spring length beyond which mangroves detach */ +#define DAMP_MANGROVE 0.1 /* damping coefficient of mangroves */ #define MANGROVE_MASS 1.5 /* mass of mangrove of radius MU */ +#define HASHX 25 /* size of hashgrid in x direction */ +#define HASHY 15 /* size of hashgrid in y direction */ +#define HASHMAX 10 /* maximal number of mangroves per hashgrid cell */ +#define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ + /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 10.0 /* max value of wave amplitude */ - -#include "hsluv.c" - #include "global_pdes.c" #include "sub_wave.c" #include "wave_common.c" double courant2, courantb2; /* Courant parameters squared */ -double circle_energy[NMAXCIRCLES]; /* energy dissipated by the circles */ -double circley_wrapped[NMAXCIRCLES]; /* position of circle centers wrapped vertically */ -double anchor_x[NMAXCIRCLES]; /* points moving circles are attached to */ -double anchor_y[NMAXCIRCLES]; /* points moving circles are attached to */ -double vx[NMAXCIRCLES]; /* x velocity of circles */ -double vy[NMAXCIRCLES]; /* y velocity of circles */ -double circlerad_initial[NMAXCIRCLES]; /* initial circle radii */ -double mass_inverse[NMAXCIRCLES]; /* inverse of mangrove mass */ -short int circle_attached[NMAXCIRCLES]; /* has value 1 if the circle is attached to its anchor */ + /*********************/ /* animation part */ @@ -411,15 +429,109 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double * evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); } +void hash_xy_to_ij(double x, double y, int ij[2]) +{ + static int first = 1; + static double lx, ly; + int i, j; + + if (first) + { + lx = XMAX - XMIN + 2.0*HASHGRID_PADDING; + ly = YMAX - YMIN + 2.0*HASHGRID_PADDING; + first = 0; + } + + i = (int)((double)HASHX*(x - XMIN + HASHGRID_PADDING)/lx); + j = (int)((double)HASHY*(y - YMIN + HASHGRID_PADDING)/ly); + + if (i<0) i = 0; + if (i>=HASHX) i = HASHX-1; + if (j<0) j = 0; + if (j>=HASHY) j = HASHY-1; + + ij[0] = i; + ij[1] = j; + +// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]); +} + + +void compute_repelling_force(int i, int j, double force[2]) +/* compute repelling force of mangrove j on mangrove i */ +{ + double x1, y1, x2, y2, distance, r, f; + + x1 = circlex[i]; + y1 = circley[i]; + x2 = circlex[j]; + y2 = circley[j]; + + distance = module2(x2 - x1, y2 - y1); + r = circlerad[i] + circlerad[j]; + if (r <= 0.0) r = 0.001*MU; + f = KREPEL/(0.001 + distance*distance); + + if ((distance > 0.0)&&(distance < REPEL_RADIUS*r)) + { + force[0] = f*(x1 - x2)/distance; + force[1] = f*(y1 - y2)/distance; + } + else + { + force[0] = 0.0; + force[1] = 0.0; + } +} + + +void update_hashgrid(int* mangrove_hashx, int* mangrove_hashy, int* hashgrid_number, int* hashgrid_mangroves) +{ + int i, j, k, n, m, ij[2], max = 0; + + printf("Updating hashgrid_number\n"); + for (i=0; i max) max = n; +// printf("Placed mangrove %i at (%i,%i) in hashgrid\n", k, ij[0], ij[1]); +// printf("%i mangroves at (%i,%i)\n", hashgrid_number[ij[0]][ij[1]], ij[0], ij[1]); + } + + printf("Maximal number of mangroves per hash cell: %i\n", max); +} + void animation() { - double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright, length; + double time, scale, diss, rgb[3], hue, y, dissip, ej, gradient[2], dx, dy, dt, xleft, xright, + length, fx, fy, force[2]; double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX]; short int *xy_in[NX], redraw = 0; - int i, j, s, ij[2]; + int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q; static int imin, imax; static short int first = 1; + + double *circleenergy, *circley_wrapped, *anchor_x, *anchor_y, *vx, *vy, *circlerad_initial, *mass_inverse; + short int *circle_attached; + int *mangrove_hashx, *mangrove_hashy, *hashgrid_number, *hashgrid_mangroves; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i= YMAX) y -= circlerad[i]; if (y <= YMIN) y += circlerad[i]; // if (y >= YMAX) y -= (YMAX - YMIN); // if (y <= YMIN) y += (YMAX - YMIN); circley_wrapped[i] = y; +// circleactive[i] = 1; if (RANDOM_RADIUS) circlerad[i] = circlerad[i]*(0.75 + 0.5*((double)rand()/RAND_MAX)); @@ -483,6 +611,9 @@ void animation() } } + /* initialise hash table for interacting mangroves */ + if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves); + if (first) /* compute box limits where circles are reset */ { /* find leftmost and rightmost circle */ @@ -515,7 +646,7 @@ void animation() for (i=0; i<=INITIAL_TIME + NSTEPS; i++) { - //printf("%d\n",i); + printf("Computing frame %d\n",i); /* compute the variance of the field to adjust color scheme */ /* the color depends on the field divided by sqrt(1 + variance) */ if (SCALE) @@ -525,83 +656,24 @@ void animation() } else scale = 1.0; + printf("Drawing wave\n"); draw_wave(phi, psi, xy_in, scale, i, PLOT); + + + printf("Evolving wave\n"); for (j=0; j 0.1*MANGROVE_EMAX) - { - dissip = 0.1*MANGROVE_EMAX; - printf("Flooring dissipation!\n"); - } - - if (circleactive[j]) - { - circle_energy[j] += dissip; - ej = circle_energy[j]; - if (ej <= MANGROVE_EMAX) - { - if (ej > 0.0) - { - hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX; - if (hue < 0.0) hue += 360.0; - } - else hue = MANGROVE_HUE_MIN; - hsl_to_rgb(hue, 0.9, 0.5, rgb); - if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue); - draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb); - - /* shrink mangrove */ - if (ej > 0.0) - { -// circlerad[j] -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX); -// if (circlerad[j] < 0.0) circlerad[j] = 0.0; - circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); - redraw = 1; - } - else circlerad[j] = circlerad_initial[j]; - } - else /* remove mangrove */ - { - circleactive[j] = 0; - /* reinitialize table xy_in */ - redraw = 1; - } - } - else /* allow disabled mangroves to recover */ - { - circle_energy[j] -= 0.15*dissip; -// circlerad[j] += 0.005*MU; -// if (circlerad[j] > MU) circlerad[j] = MU; -// if ((circle_energy[j] < 0.0)&&(circlerad[j] > 0.0)) - if (circle_energy[j] < 0.0) - { - circleactive[j] = 1; -// circlerad[j] = circlerad[j]*(0.75 + 0.5*((double)rand()/RAND_MAX)); - circlerad[j] = circlerad_initial[j]; - circle_energy[j] = -MANGROVE_EMAX; - /* reinitialize table xy_in */ - redraw = 1; - } - - } - -// printf("Circle %i, energy %.5lg\n", j, circle_energy[j]); - } /* move mangroves */ if (MOVE_MANGROVES) for (j=0; j ", j, circlex[j], circley[j]); @@ -617,11 +689,41 @@ void animation() dy += DT_MANGROVE*(-KSPRING*(circley_wrapped[j] - anchor_y[j])); } + /* compute repelling force from other mangroves */ + if (REPELL_MANGROVES) + { + /* determine neighboring grid points */ + i0 = mangrove_hashx[j]; + iminus = i0 - 1; if (iminus < 0) iminus = 0; + iplus = i0 + 1; if (iplus >= HASHX) iplus = HASHX-1; + + j0 = mangrove_hashy[j]; + jminus = j0 - 1; if (jminus < 0) jminus = 0; + jplus = j0 + 1; if (jplus >= HASHY) jplus = HASHY-1; + + fx = 0.0; + fy = 0.0; + for (p=iminus; p<= iplus; p++) + for (q=jminus; q<= jplus; q++) + for (k=0; k 0.001) printf("Force on mangrove %i: (%.3f, %.3f)\n", j, fx, fy); + + dx += DT_MANGROVE*fx; + dy += DT_MANGROVE*fy; + } + /* detach mangrove if spring is too long */ if (DETACH_MANGROVES) { length = module2(circlex[j] - anchor_x[j], circley_wrapped[j] - anchor_y[j]); - if (j%NGRIDY == 0) printf("spring length %.i: %.3lg\n", j, length); +// if (j%NGRIDY == 0) printf("spring length %.i: %.3lg\n", j, length); // if (length > L_DETACH) circle_attached[j] = 0; if (length*mass_inverse[j] > L_DETACH) circle_attached[j] = 0; } @@ -638,9 +740,9 @@ void animation() circlex[j] += vx[j]*DT_MANGROVE; circley[j] += vy[j]*DT_MANGROVE; circley_wrapped[j] += vy[j]*DT_MANGROVE; - if (j%NGRIDY == 0) - printf("circle %.i: (dx,dy) = (%.3lg,%.3lg), (vx,vy) = (%.3lg,%.3lg)\n", - j, circlex[j]-anchor_x[j], circley[j]-anchor_y[j], vx[j], vy[j]); +// if (j%NGRIDY == 0) +// printf("circle %.i: (dx,dy) = (%.3lg,%.3lg), (vx,vy) = (%.3lg,%.3lg)\n", +// j, circlex[j]-anchor_x[j], circley[j]-anchor_y[j], vx[j], vy[j]); } else { @@ -655,10 +757,148 @@ void animation() if (circley_wrapped[j] >= YMAX) circley_wrapped[j] = YMAX; // if (j%NGRIDY == 0) printf("(%.3lg, %.3lg)\n", circlex[j], circley[j]); - + redraw = 1; } + /* test for debugging */ + if (1) for (j=0; j 0.1*MANGROVE_EMAX) + { + dissip = 0.1*MANGROVE_EMAX; + printf("Flooring dissipation!\n"); + } + + if (circleactive[j]) + { + circleenergy[j] += dissip; + ej = circleenergy[j]; +// printf("ej = %.3f\n", ej); + if (ej <= MANGROVE_EMAX) + { + if (ej > 0.0) + { + hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX; + if (hue < 0.0) hue += 360.0; + } + else hue = MANGROVE_HUE_MIN; + hsl_to_rgb(hue, 0.9, 0.5, rgb); +// if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue); + draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb); + + /* shrink mangrove */ + if ((ERODE_MANGROVES)&&(ej > 0.0)) + { + circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); + redraw = 1; + } + else circlerad[j] = circlerad_initial[j]; + + } + else /* remove mangrove */ + { + circleactive[j] = 0; + /* reinitialize table xy_in */ + redraw = 1; + } + } + else if (RECOVER_MANGROVES) /* allow disabled mangroves to recover */ + { + circleenergy[j] -= 0.15*dissip; + printf("Circle %i energy %.3lg\n", j, circleenergy[j]); + if (circleenergy[j] < 0.0) + { + printf("Reactivating circle %i?\n", j); + /* THE PROBLEM occurs when circleactive[0] is set to 1 again */ + if (j>0) circleactive[j] = 1; + circlerad[j] = circlerad_initial[j]; + circleenergy[j] = -MANGROVE_EMAX; + /* reinitialize table xy_in */ + redraw = 1; + } + + } + } + + /* compute energy dissipated in obstacles */ +/* if (ERODE_MANGROVES) for (j=0; j 0.1*MANGROVE_EMAX) +// { +// dissip = 0.1*MANGROVE_EMAX; +// printf("Flooring dissipation!\n"); +// } +// +// if (circleactive[j]) +// { +// circleenergy[j] += dissip; +// ej = circleenergy[j]; +// printf("ej = %.3f\n", ej); +// if (ej <= MANGROVE_EMAX) +// { +// if (ej > 0.0) +// { +// hue = MANGROVE_HUE_MIN + (MANGROVE_HUE_MAX - MANGROVE_HUE_MIN)*ej/MANGROVE_EMAX; +// if (hue < 0.0) hue += 360.0; +// } +// else hue = MANGROVE_HUE_MIN; +// hsl_to_rgb(hue, 0.9, 0.5, rgb); +// if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue); +// draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb); +// +// /* shrink mangrove */ +// if (ej > 0.0) +// { +// circlerad[j] -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX); +// if (circlerad[j] < 0.0) circlerad[j] = 0.0; +// circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); +// redraw = 1; +// } +// else circlerad[j] = circlerad_initial[j]; +// } +// else /* remove mangrove */ +// { +// circleactive[j] = 0; + /* reinitialize table xy_in */ +// redraw = 1; +// } +// } +// else /* allow disabled mangroves to recover */ +// { +// circleenergy[j] -= 0.15*dissip; +// printf("ej = %.3f\n", circleenergy[j]); +// circlerad[j] += 0.005*MU; +// if (circlerad[j] > MU) circlerad[j] = MU; +// if ((circleenergy[j] < 0.0)&&(circlerad[j] > 0.0)) +// if (circleenergy[j] < 0.0) +// { +// circleactive[j] = 1; +// circlerad[j] = circlerad[j]*(0.75 + 0.5*((double)rand()/RAND_MAX)); +// circlerad[j] = circlerad_initial[j]; +// circleenergy[j] = -MANGROVE_EMAX; + /* reinitialize table xy_in */ +// redraw = 1; +// } + +// } + +// printf("Circle %i, energy %.5lg\n", j, circleenergy[j]); +// } + + printf("Updating hashgrid\n"); + if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves); + + + printf("Drawing billiard\n"); draw_billiard(); glutSwapBuffers(); @@ -666,7 +906,7 @@ void animation() if (redraw) { printf("Reinitializing xy_in\n"); - init_xyin_xrange(xy_in, imin, NX-1); + init_xyin_xrange(xy_in, imin, NX); // init_xyin_xrange(xy_in, imin, imax); } redraw = 0; @@ -701,6 +941,21 @@ void animation() free(psi_tmp[i]); free(xy_in[i]); } + + free(circleenergy); + free(circley_wrapped); + free(anchor_x); + free(anchor_y); + free(vx); + free(vy); + free(circlerad_initial); + free(mass_inverse); + free(circle_attached); + + free(mangrove_hashx); + free(mangrove_hashy); + free(hashgrid_number); + free(hashgrid_mangroves); } diff --git a/particle_billiard.c b/particle_billiard.c index d56e26b..5a183c3 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -32,13 +32,18 @@ #define MOVIE 0 /* set to 1 to generate movie */ -#define WINWIDTH 1280 /* window width */ +// #define WINWIDTH 1280 /* window width */ +#define WINWIDTH 720 /* window width */ #define WINHEIGHT 720 /* window height */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#define XMIN -1.4 +#define XMAX 1.4 /* x interval */ +#define YMIN -1.4 +#define YMAX 1.4 /* y interval for 9/16 aspect ratio */ +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ @@ -60,7 +65,7 @@ // #define LAMBDA 1.4 /* parameter controlling shape of domain */ // #define MU 0.2 /* second parameter controlling shape of billiard */ -#define LAMBDA 1.5 /* parameter controlling shape of domain */ +#define LAMBDA 1.3 /* parameter controlling shape of domain */ #define MU 0.3 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NPOLY 4 /* number of sides of polygon */ @@ -74,15 +79,16 @@ /* Simulation parameters */ -#define NPART 10000 /* number of particles */ +#define NPART 100 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ #define CYCLE 1 /* set to 1 for closed curve (start in all directions) */ -#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ +#define SHOWTRAILS 1 /* set to 1 to keep trails of the particles */ +#define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */ -#define NSTEPS 3500 /* number of frames of movie */ -#define TIME 1200 /* time between movie frames, for fluidity of real-time simulation */ +#define NSTEPS 1300 /* number of frames of movie */ +#define TIME 2500 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.00001 /* integration step */ #define NVID 150 /* number of iterations between images displayed on screen */ @@ -94,12 +100,14 @@ /* Colors and other graphical parameters */ -#define NCOLORS 16 /* number of colors */ +#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ + +#define NCOLORS 64 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ -#define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ +#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ -#define LENGTH 0.01 /* length of velocity vectors */ +#define LENGTH 0.03 /* length of velocity vectors */ #define BILLIARD_WIDTH 2 /* width of billiard */ #define PARTICLE_WIDTH 2 /* width of particles */ #define FRONT_WIDTH 3 /* width of wave front */ @@ -170,6 +178,31 @@ void init_drop_config(double x0, double y0, double angle1, double angle2, double } } +void init_partial_drop_config(double x0, double y0, double angle1, double angle2, int particle1, int particle2, + int col, double *configs[NPARTMAX], int color[NPARTMAX], int newcolor[NPARTMAX]) +/* initialize configuration: drop at (x0,y0) for a range of particles */ +{ + int i; + double dalpha, alpha; + double conf[2], pos[2]; + + while (angle2 < angle1) angle2 += DPI; + if (particle2 - particle1 > 1) dalpha = (angle2 - angle1)/((double)(particle2 - particle1-1)); + else dalpha = 0.0; + for (i=particle1; i= NCOLORS) color[i] -= NCOLORS; +// } +// } + + configs[i][2] += DPHI; + + cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; + sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; + len = configs[i][2] + LENGTH; + if (len > configs[i][3]) len = configs[i][3]; + + x0 = configs[i][4]; + y0 = configs[i][5]; + x1 = configs[i][4] + configs[i][2]*cosphi; + y1 = configs[i][5] + configs[i][2]*sinphi; + x2 = configs[i][4] + len*cosphi; + y2 = configs[i][5] + len*sinphi; + + /* test whether particle does not escape billiard */ + if ((TEST_ACTIVE)&&(active[i])) active[i] = xy_in_billiard(x1, y1); + + if (active[i]) + { + rgb_color_scheme(color[i], rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + glBegin(GL_LINE_STRIP); + glVertex2d(SCALING_FACTOR*x0, SCALING_FACTOR*y0); + glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*y2); + glEnd (); + } + +// if (configs[i][2] > configs[i][3] - DPHI) +// { +// glBegin(GL_LINE_STRIP); +// glVertex2d(SCALING_FACTOR*x0, SCALING_FACTOR*y0); +// glVertex2d(SCALING_FACTOR*configs[i][6], SCALING_FACTOR*configs[i][7]); +// glEnd (); +// } + } + if (DRAW_BILLIARD) draw_billiard(); +} + + void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) /* draw the particles */ { @@ -332,7 +430,8 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a { // printf("reflecting particle %i\n", i); c = vbilliard(configs[i]); - if (c>=0) color[i]++; +// if (c>=0) color[i]++; + if ((!RAINBOW_COLOR)&&(c>=0)) color[i]++; if (!RAINBOW_COLOR) { color[i]++; @@ -354,7 +453,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a void animation() { - double time, dt, alpha, r; + double time, dt, alpha, r, rgb[3]; double *configs[NPARTMAX]; int i, j, resamp = 1, s, i1, i2; int *color, *newcolor, *active; @@ -430,15 +529,27 @@ void animation() } sleep(SLEEP1); + + + /* initialize drops in different colors */ + init_partial_drop_config(0.0, 0.0, 0.0, DPI, 0, 2*NPART/5, 0, configs, color, newcolor); + init_partial_drop_config(0.0, 0.8, 0.0, DPI, 2*NPART/5, 4*NPART/5, 10, configs, color, newcolor); + init_partial_drop_config(1.2, 0.1, 0.0, DPI, 4*NPART/5, NPART, 36, configs, color, newcolor); for (i=0; i<=NSTEPS; i++) { graph_movie(TIME, newcolor, configs, active); - draw_config(newcolor, configs, active); + if (SHOWTRAILS) draw_config_showtrails(newcolor, configs, active); + else draw_config(newcolor, configs, active); +// draw_config(newcolor, configs, active); if (DRAW_BILLIARD) draw_billiard(); for (j=0; j= 0) ncol++; + if ((c >= 0)&&(circlecolor[c] == 0)) nobst++; circlecolor[c]++; + + /* take care of circles doubled because of periodic boundary conditions */ + if ((circleactive[c])&&(B_DOMAIN == D_CIRCLES_IN_TORUS)&&(partner_circle[c] != c)) circlecolor[partner_circle[c]]++; + newcircle[c] = 10; /* update free path statistics */ - k = (int)((double)NPATHBINS*configs[i][3]/PATHLMAX); - if (k < NPATHBINS) npath[k]++; - -// if (circlecolor[c] > nmaxpeg) nmaxpeg = circlecolor[c]; - - -// if (circlecolor[c] > NCOLORS) circlecolor[c] = 1; -// if (c>=0) color[i]++; - + if (ncol > 1) /* disregard very first collision */ + { + if (B_DOMAIN != D_CIRCLES_IN_TORUS) + { + k = (int)((double)NPATHBINS*configs[i][3]/PATHLMAX); + if (k < NPATHBINS) npath[k]++; + if (total_pathlength > max_free_path) max_free_path = total_pathlength; + total_pathlength = 0.0; + } + else /* case with periodic boundary conditions */ + { + if (c >= 0) /* a circle is hit, update histogram */ + { +// printf("total path length %.3lg\n", total_pathlength); + k = (int)((double)NPATHBINS*total_pathlength/PATHLMAX); + if (k < NPATHBINS) npath[k]++; + if (total_pathlength > max_free_path) max_free_path = total_pathlength; + total_pathlength = 0.0; + } + } + } if (!RAINBOW_COLOR) { @@ -514,8 +558,10 @@ void print_particle_numbers(double *configs[NPARTMAX]) void draw_statistics() { - int i, n, colmax = 35, pegcollisions[70], nypegs = 70, meanpegs = 0, total_coll = 0, ymax = 0, meanbins = 0, total_bin = 0; - double x, y, yscale = 120.0, y0, dx, rgb[3], xshift; + int i, n, colmax = 55, pegcollisions[90], nypegs = 70, meanpegs = 0, meansquarepegs = 0, total_coll = 0, ymax = 0, + meanbins = 0, meansquarebins, total_bin = 0, stephits = 10, tickstephits = 1; + double x, y, yscale = 110.0, y0, dx, rgb[3], xshift, coll_mean, coll_stdv, path_mean, path_stdv, ebin, len_over_bin, + steppath = 0.5; char message[50]; glLineWidth(1); @@ -525,10 +571,10 @@ void draw_statistics() xshift = XMIN + 0.3; rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 1.0; - /* histogramm of number of collisions per peg */ + /* histogram of number of collisions per peg */ for (i=0; i ymax) ymax = npath[i]; total_bin += npath[i]; meanbins += i*npath[i]; + meansquarebins += i*i*npath[i]; } yscale = 0.9*(YMAX-y0)*(double)ymax; @@ -620,13 +681,13 @@ void draw_statistics() x = xshift + (double)i*dx; y = y0 + (double)npath[i]*YMAX/yscale; - erase_rectangle(x, y0, x+dx, y, rgb); + if (y > y0) erase_rectangle(x, y0, x+dx, y, rgb); } glColor3f(1.0, 1.0, 1.0); glBegin(GL_LINE_STRIP); - /* histogramm */ - for (i=1; i 0) { x = xshift + (double)i*dx; y = y0 + (double)npath[i]*YMAX/yscale; @@ -636,19 +697,23 @@ void draw_statistics() glVertex2d(x+dx, y0); glVertex2d(x, y0); } - glVertex2d(xshift, y0); + glEnd (); + + glBegin(GL_LINE_STRIP); glVertex2d(xshift, YMAX - 0.1); + glVertex2d(xshift, y0); + glVertex2d(xshift + (double)NPATHBINS*dx, y0); glEnd (); - for (x = 0.5; x < PATHLMAX; x+=0.5) + for (x = steppath; x < PATHLMAX; x+=steppath) { i = (int)(x*(double)NPATHBINS/PATHLMAX); sprintf(message, "%.2f", x); write_text_fixedwidth(xshift + (double)i*dx - 0.1, y0 - 0.12, message); } - for (x = 0.1; x < PATHLMAX; x+=0.1) + for (x = steppath; x < PATHLMAX; x+=steppath) { i = (int)(x*(double)NPATHBINS/PATHLMAX); glBegin(GL_LINE_STRIP); @@ -656,12 +721,24 @@ void draw_statistics() glVertex2d(xshift + (double)i*dx, y0 + 0.025); glEnd (); } - + + + ebin = (double)meanbins/(double)total_bin; /* mean bin */ + len_over_bin = PATHLMAX/(double)NPATHBINS; /* conversion from bin to path length */ + path_mean = ebin*len_over_bin; /* mean free path */ + path_stdv = sqrt((double)meansquarebins/(double)total_bin - ebin*ebin)*len_over_bin; + sprintf(message, "free path"); write_text_fixedwidth(XMAX - 0.6, y0 - 0.12, message); - sprintf(message, "Mean free path %.4lg", (double)meanbins*PATHLMAX/(double)(total_bin*NPATHBINS)); + sprintf(message, "Max free path %.4lg", max_free_path); write_text(2.2, YMAX - 0.3, message); + + sprintf(message, "Mean free path %.4lg", path_mean); + write_text(2.2, YMAX - 0.5, message); + + sprintf(message, "Stdv free path %.4lg", path_stdv); + write_text(2.2, YMAX - 0.7, message); } @@ -681,24 +758,38 @@ void animation(int circle_config) configs[i] = (double *)malloc(8*sizeof(double)); /* init circle configuration if the domain is D_CIRCLES */ - if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN)) + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN) + ||(B_DOMAIN == D_CIRCLES_IN_TORUS)) init_circle_config_pinball(circle_config); /* remove discs that are not in domain */ - if ((B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN)) + if ((B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN)) +// ||(B_DOMAIN == D_CIRCLES_IN_TORUS)) for (i=0; i 0.99) circleactive[i] = 0; if (vabs(circlex[i]) + circlerad[i] > 0.99*LAMBDA) circleactive[i] = 0; } + else if (B_DOMAIN == D_CIRCLES_IN_TORUS) + for (i=0; i 1.0) circleactive[i] = 0; + if (vabs(circlex[i]) - circlerad[i] > LAMBDA) circleactive[i] = 0; + } + +// for (i=0; i 1.0) circleactive[i] = 0; /* initialize system by putting particles in a given point with a range of velocities */ r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY); +// init_drop_config(0.4, 0.0, PID, 0.5*PID, configs); + init_drop_config(0.0, 0.0, -0.45*PID, 0.45*PID, configs); + // init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); - init_drop_config(0.5, 0.1, -0.5*PID, 0.5*PID, configs); -// init_drop_config(0.0, 0.0, -0.5*PID, 0.5*PID, configs); +// init_drop_config(0.5, 0.1, -0.5*PID, 0.5*PID, configs); // init_drop_config(-1.4, 0.0, -0.5*PID, 0.5*PID, configs); // init_drop_config(0.5, 0.5, -1.0, 1.0, configs); // init_sym_drop_config(-1.0, 0.5, -PID, PID, configs); @@ -795,9 +886,9 @@ void animation(int circle_config) for (j=0; j nmaxpeg) nmaxpeg = circlecolor[j]; - sprintf(message, "max hits per peg: %d", nmaxpeg); - glColor3f(1.0, 1.0, 1.0); - write_text(-0.6, YMIN + 0.08, message); +// sprintf(message, "max hits per peg: %d", nmaxpeg); +// glColor3f(1.0, 1.0, 1.0); +// write_text(-0.6, YMIN + 0.08, message); // write_text(-0.3, YMIN + 0.04, message); draw_statistics(); @@ -849,8 +940,8 @@ void display(void) } animation(CIRCLE_PATTERN); - // animation(C_TRI); + // animation(C_GOLDEN_SPIRAL); // animation(C_SQUARE); // animation(C_HEX); diff --git a/schrodinger.c b/schrodinger.c index 6cb47b5..c4cdecc 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -39,10 +39,12 @@ /* General geometrical parameters */ -#define WINWIDTH 1280 /* window width */ +// #define WINWIDTH 1280 /* window width */ +#define WINWIDTH 720 /* window width */ #define WINHEIGHT 720 /* window height */ -#define NX 1280 /* number of grid points on x axis */ +// #define NX 1280 /* number of grid points on x axis */ +#define NX 720 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ // #define NX 640 /* number of grid points on x axis */ // #define NY 360 /* number of grid points on y axis */ @@ -52,22 +54,24 @@ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#define YMIN -2.0 +#define YMAX 2.0 /* y interval for 9/16 aspect ratio */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table, see list in global_pdes.c */ -#define B_DOMAIN 9 /* choice of domain shape */ +#define B_DOMAIN 19 /* choice of domain shape */ #define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ -#define MU 0.05 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */ +#define MU 1.75 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 3 /* depth of computation of Menger gasket */ @@ -78,13 +82,25 @@ #define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard in sub_wave.c */ /* Physical patameters of wave equation */ // #define DT 0.00000005 -#define DT 0.000000005 +#define DT 0.00000001 +// #define DT 0.000000005 // #define DT 0.000000005 #define HBAR 1.0 @@ -94,8 +110,9 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 200 /* number of frames of movie */ -#define NVID 1200 /* number of iterations between images displayed on screen */ +#define NSTEPS 1400 /* number of frames of movie */ +#define NVID 2000 /* number of iterations between images displayed on screen */ +// #define NVID 1200 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ @@ -111,7 +128,7 @@ /* Plot type, see list in global_pdes.c */ -#define PLOT 10 +#define PLOT 11 /* Color schemes, see list in global_pdes.c */ @@ -133,8 +150,6 @@ #define HUEMEAN 150.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -150.0 /* amplitude of variation of hue for color scheme C_HUE */ -#include "hsluv.c" - #include "global_pdes.c" #include "sub_wave.c" @@ -401,7 +416,7 @@ void animation() printf("Integration step %.3lg\n", intstep); /* initialize wave wave function */ - init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.25, phi, psi, xy_in); + init_coherent_state(0.5, 0.0, 40.0, 0.0, 0.25, phi, psi, xy_in); // init_coherent_state(0.0, 0.0, 0.0, 5.0, 0.03, phi, psi, xy_in); // init_coherent_state(-0.5, 0.0, 1.0, 1.0, 0.05, phi, psi, xy_in); diff --git a/sub_part_billiard.c b/sub_part_billiard.c index d3c36b6..e83df9a 100644 --- a/sub_part_billiard.c +++ b/sub_part_billiard.c @@ -1,5 +1,8 @@ +#include "colormaps.c" + #define DUMMY_ABSORBING -1000.0 /* dummy value of config[0] for absorbing circles */ #define BOUNDARY_SHIFT 100000.0 /* shift of boundary parametrisation for circles in domain */ +#define DUMMY_SIDE_ABS -10000 /* dummy value of returned side for absorbing circles */ long int global_time = 0; /* counter to keep track of global time of simulation */ int nparticles=NPART; @@ -133,44 +136,6 @@ void init() /* initialisation of window */ } -void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */ -/* h = hue, s = saturation, l = luminosity */ -{ - double c = 0.0, m = 0.0, x = 0.0; - - c = (1.0 - fabs(2.0 * l - 1.0)) * s; - m = 1.0 * (l - 0.5 * c); - x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0)); - - if (h >= 0.0 && h < 60.0) - { - rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m; - } - else if (h < 120.0) - { - rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m; - } - else if (h < 180.0) - { - rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m; - } - else if (h < 240.0) - { - rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m; - } - else if (h < 300.0) - { - rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m; - } - else if (h < 360.0) - { - rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m; - } - else - { - rgb[0] = m; rgb[1] = m; rgb[2] = m; - } -} void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ { @@ -186,6 +151,7 @@ void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ /* saturation = r, luminosity = 0.5 */ } + void rgb_color_scheme_lum(int i, double lum, double rgb[3]) /* color scheme */ { double hue, y, r; @@ -349,6 +315,18 @@ void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) } +void draw_initial_condition_circle(double x, double y, double r, int color) +/* draws a colored circle to mark initial condition */ +{ + double rgb[3]; + + rgb_color_scheme(color, rgb); + draw_colored_circle(x, y, r, NSEG, rgb); + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + glColor3f(rgb[0], rgb[1], rgb[2]); + glLineWidth(4); + draw_circle(x, y, r, NSEG); +} int in_circle(double x, double y, double r) /* test whether (x,y) is in circle of radius r */ @@ -1163,10 +1141,13 @@ void draw_billiard() /* draws the billiard boundary */ { if (circleactive[k] == 2) { - hsl_to_rgb(150.0, 0.9, 0.4, rgb); - glColor3f(rgb[0], rgb[1], rgb[2]); +// hsl_to_rgb(150.0, 0.9, 0.4, rgb); +// glColor3f(rgb[0], rgb[1], rgb[2]); + glColor3f(0.0, 1.0, 0.0); + rgb[0] = 0.0; rgb[1] = 0.9; rgb[2] = 0.0; + draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb); } - draw_circle(circlex[k], circley[k], circlerad[k], NSEG); + else draw_circle(circlex[k], circley[k], circlerad[k], NSEG); init_billiard_color(); } @@ -1227,6 +1208,55 @@ void draw_billiard() /* draws the billiard boundary */ glEnd (); break; } + case (D_CIRCLES_IN_TORUS): + { + rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; + for (k=0; k= 1) + { + rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); + newcircle[k]--; + } + else rgb_color_scheme(circlecolor[k], rgb1); + draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); + } + } + init_billiard_color(); + for (k=0; k= 1) + { + if (circleactive[k] == 2) + { + hsl_to_rgb(150.0, 0.9, 0.4, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + } + draw_circle(circlex[k], circley[k], circlerad[k], NSEG); + init_billiard_color(); + } + + /* draw shooter position for laser pattern */ + if (CIRCLE_PATTERN == C_LASER) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circlerad[ncircles-1], NSEG); + } + + init_billiard_color(); + + glBegin(GL_LINE_LOOP); + glVertex2d(LAMBDA, -1.0); + glVertex2d(LAMBDA, 1.0); + glVertex2d(-LAMBDA, 1.0); + glVertex2d(-LAMBDA, -1.0); + glEnd(); + break; + } default: { printf("Function draw_billiard not defined for this billiard \n"); @@ -3924,7 +3954,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ int vcircles_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ { - double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10; + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1000.0; double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; int i, nt = 0, nscat[ncircles], ntmin; @@ -3991,8 +4021,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ } else /* there is no intersection - set dummy values */ { - config[0] = 0.0; - config[1] = 0.0; + config[0] = DUMMY_ABSORBING; + config[1] = PI; config[2] = 0.0; config[3] = rlarge; config[4] = pos[0]; /* start position */ @@ -4012,7 +4042,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ c = pos_circles(config, pos, &alpha); - vcircles_xy(config, alpha, pos); + c = vcircles_xy(config, alpha, pos); return(c); } @@ -4125,9 +4155,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { config[0] = DUMMY_ABSORBING; config[1] = PI; +// return(DUMMY_SIDE_ABS); } - return(nscat[ntmin]); + if (ABSORBING_CIRCLES) return(DUMMY_SIDE_ABS); + else return(nscat[ntmin]); } else /* there is no intersection with the circles - compute intersection with boundary */ { @@ -4150,13 +4182,14 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ c = pos_circles_in_rect(config, pos, &alpha); - vcircles_in_rect_xy(config, alpha, pos); +// vcircles_in_rect_xy(config, alpha, pos); + c = vcircles_in_rect_xy(config, alpha, pos); return(c); } - /****************************************************************************************/ +/****************************************************************************************/ /* billiard with circular scatterers in a genus n surface (polygon with identifies opposite sides) */ /****************************************************************************************/ @@ -4298,6 +4331,151 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ } +/****************************************************************************************/ +/* billiard with circular scatterers in a torus (rectangle with periodic boundary conditions) */ +/****************************************************************************************/ + + int pos_circles_in_torus(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of circle */ + /* position varies between 0 and ncircles*2Pi for circles and between -BOUNDARY_SHIFT and 0 for boundary*/ + /* returns number of hit circle */ + { + double angle; + int ncirc, c; + + if (conf[0] >= 0) + { + ncirc = (int)(conf[0]/DPI); + if (ncirc >= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle); + pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + else /* particle starts on boundary */ + { +// conf[0] += 4.0*(LAMBDA + 1.0); + conf[0] += BOUNDARY_SHIFT; + c = pos_rectangle(conf, pos, alpha); + +// conf[0] -= 4.0*(LAMBDA + 1.0); + conf[0] -= BOUNDARY_SHIFT; + + return(-c-1); + } + } + + + int vcircles_in_torus_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ +// double config[8], alpha, pos[2]; + + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, nt = 0, nscat[ncircles], ntmin, side; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circlex[i], yint[nt] - circley[i]); + + /* test wether intersection is in rectangle */ + if ((vabs(xint[nt]) < LAMBDA)&&(vabs(yint[nt]) < 1.0)) nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + + return(nscat[ntmin]); + } + else /* there is no intersection with the circles - compute intersection with boundary */ + { + + side = vrectangle_xy(config, alpha, pos); + + if (config[0] < 2.0*LAMBDA) config[0] = 4.0*LAMBDA + 2.0 - config[0]; + else if (config[0] < 2.0*LAMBDA + 2.0) config[0] = 6.0*LAMBDA + 4.0 - config[0]; + else if (config[0] < 4.0*LAMBDA + 2.0) config[0] = 4.0*LAMBDA + 2.0 - config[0]; + else config[0] = 6.0*LAMBDA + 4.0 - config[0]; + + config[0] -= BOUNDARY_SHIFT; + config[1] = PI - config[1]; +// config[0] -= 4.0*(LAMBDA+1.0); + +// printf("Hit side %i\n", side); +// print_config(config); + return(side - 5); + } + } + + int vcircles_in_torus(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles_in_torus(config, pos, &alpha); + + vcircles_in_torus_xy(config, alpha, pos); + + return(c); + } + /****************************************************************************************/ /* general billiard */ @@ -4397,6 +4575,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(pos_circles_in_genusn(conf, pos, alpha)); break; } + case (D_CIRCLES_IN_TORUS): + { + return(pos_circles_in_torus(conf, pos, alpha)); + break; + } default: { printf("Function pos_billiard not defined for this billiard \n"); @@ -4500,6 +4683,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vcircles_in_genusn_xy(config, alpha, pos)); break; } + case (D_CIRCLES_IN_TORUS): + { + return(vcircles_in_torus_xy(config, alpha, pos)); + break; + } default: { printf("Function vbilliard_xy not defined for this billiard \n"); @@ -4642,6 +4830,13 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vcircles_in_genusn(config)); break; } + case (D_CIRCLES_IN_TORUS): + { + c = pos_circles_in_torus(config, pos, &alpha); + + return(vcircles_in_torus(config)); + break; + } default: { printf("Function vbilliard not defined for this billiard \n"); @@ -4843,6 +5038,18 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ } break; } + case D_CIRCLES_IN_TORUS: /* same as D_CIRCLES_IN_RECT */ + { + if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0); + else + { + condition = 1; + for (k=0; k LAMBDA) xx[1] = 2.0*LAMBDA - xx[1]; + xx[2] = -xx[0]; + xx[3] = -xx[1]; + + yy[0] = 0.5*(y_shooter + y_target); + yy[1] = 1.0 - 0.5*(y_target - y_shooter); + if (yy[1] > 1.0) yy[1] = 2.0 - yy[1]; + yy[2] = -yy[0]; + yy[3] = -yy[1]; + +// xx[0] = 0.5*(x_shooter + x_target); +// xx[1] = LAMBDA - 0.5*(x_target - x_shooter); // xx[2] = -xx[0]; // xx[3] = -xx[1]; // -// yy[0] = 0.5*(Y_SHOOTER + Y_TARGET); -// yy[1] = 1.0 - 0.5*(Y_TARGET - Y_SHOOTER); +// yy[0] = 0.5*(y_shooter + y_target); +// yy[1] = 1.0 - 0.5*(y_target - y_shooter); // yy[2] = -yy[0]; // yy[3] = -yy[1]; -// -// for (i = 0; i < 4; i++) -// for (j = 0; j < 4; j++) -// { -// circlex[4*i + j] = xx[i]; -// circley[4*i + j] = yy[j]; -// -// } -// -// circlex[ncircles - 1] = X_TARGET; -// circley[ncircles - 1] = Y_TARGET; -// -// for (i=0; i YMAX) y -= height; + circlerad[n] = MU; + circleactive[n] = 1; + circlecolor[n] = 0; + newcircle[n] = 0; + } + break; + } + case (C_GOLDEN_SPIRAL): + { + ncircles = 1; + circlex[0] = 0.0; + circley[0] = 0.0; + + gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ + phi = 0.0; + r0 = 2.0*MU; + r = r0 + MU; + + for (i=0; i YMIN - MU)) circleactive[i] = 1; + } + break; + } + case (C_RAND_DISPLACED): + { + ncircles = NCX*NCY; + dy = (YMAX - YMIN)/((double)NCY); + for (i = 0; i < NCX; i++) + for (j = 0; j < NCY; j++) + { + n = NCY*i + j; + circlex[n] = ((double)(i-NCX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circley[n] = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circlerad[n] = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circleactive[n] = 1; + circlecolor[n] = 0; + newcircle[n] = 0; + } + break; + } + case (C_RAND_POISSON): + { + ncircles = NPOISSON; + for (n = 0; n < NPOISSON; n++) + { + circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circley[n] = (BOXYMAX - BOXYMIN)*(double)rand()/RAND_MAX + BOXYMIN; + circlerad[n] = MU; + circleactive[n] = 1; + circlecolor[n] = 0; + newcircle[n] = 0; + double_circle[n] = 0; + partner_circle[n] = n; + + /* take care of periodic boundary conditions */ + if (B_DOMAIN == D_CIRCLES_IN_TORUS) + { + /* inactivate circles in corners */ + if ((vabs(circlex[n]) > LAMBDA - MU)&&((circley[n] < BOXYMIN + MU)||(circley[n] > BOXYMAX - MU))) + circleactive[n] = 0; + + if (circlex[n] < - LAMBDA + MU) + { + circlex[ncircles] = circlex[n] + 2.0*LAMBDA; + circley[ncircles] = circley[n]; + partner_circle[ncircles] = n; + partner_circle[n] = ncircles; + ncircles++; + } + else if (circlex[n] > LAMBDA - MU) + { + circlex[ncircles] = circlex[n] - 2.0*LAMBDA; + circley[ncircles] = circley[n]; + partner_circle[ncircles] = n; + partner_circle[n] = ncircles; + ncircles++; + } + if (circley[n] < BOXYMIN + MU) + { + circlex[ncircles] = circlex[n]; + circley[ncircles] = circley[n] + BOXYMAX - BOXYMIN; + partner_circle[ncircles] = n; + partner_circle[n] = ncircles; + ncircles++; + } + else if (circley[n] > BOXYMAX - MU) + { + circlex[ncircles] = circlex[n]; + circley[ncircles] = circley[n] - BOXYMAX + BOXYMIN; + partner_circle[ncircles] = n; + partner_circle[n] = ncircles; + ncircles++; + } + } + } + printf("%i circles\n", ncircles); + if (B_DOMAIN == D_CIRCLES_IN_TORUS) for (n = NPOISSON; n < ncircles; n++) + { +// printf("circle %i at (%.3f, %.3f)\n", n, circlex[n], circley[n]); + circlerad[n] = MU; + if (circleactive[partner_circle[n]]) circleactive[n] = 1; + circlecolor[n] = 0; + newcircle[n] = 0; + double_circle[n] = 1; + } + + break; + } + case (C_POISSON_DISC): + { + printf("Generating Poisson disc sample\n"); + /* generate first circle */ + circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + active_poisson[0] = 1; + n_p_active = 1; + ncircles = 1; + + while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) + { + /* randomly select an active circle */ + i = rand()%(ncircles); + while (!active_poisson[i]) i = rand()%(ncircles); +// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN); + } + if (far) /* accept new circle */ + { + printf("New circle at (%.3f,%.3f) accepted\n", x, y); + circlex[ncircles] = x; + circley[ncircles] = y; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + circlecolor[ncircles] = 0; + newcircle[ncircles] = 0; + active_poisson[ncircles] = 1; + ncircles++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { +// printf("No candidates work, inactivate circle %i\n", i); + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + break; + } +// case (C_LASER): +// { +// ncircles = 17; +// +// xx[0] = 0.5*(X_SHOOTER + X_TARGET); +// xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER); +// xx[2] = -xx[0]; +// xx[3] = -xx[1]; +// +// yy[0] = 0.5*(Y_SHOOTER + Y_TARGET); +// yy[1] = 1.0 - 0.5*(Y_TARGET - Y_SHOOTER); +// yy[2] = -yy[0]; +// yy[3] = -yy[1]; +// +// for (i = 0; i < 4; i++) +// for (j = 0; j < 4; j++) +// { +// circlex[4*i + j] = xx[i]; +// circley[4*i + j] = yy[j]; +// +// } +// +// circlex[ncircles - 1] = X_TARGET; +// circley[ncircles - 1] = Y_TARGET; +// +// for (i=0; i= 0.0 && h < 60.0) - { - rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m; - } - else if (h < 120.0) - { - rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m; - } - else if (h < 180.0) - { - rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m; - } - else if (h < 240.0) - { - rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m; - } - else if (h < 300.0) - { - rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m; - } - else if (h < 360.0) - { - rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m; - } - else - { - rgb[0] = m; rgb[1] = m; rgb[2] = m; - } -} - -void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* color conversion from HSL to RGB */ -{ - double r, g, b; - - switch (COLOR_PALETTE) { - case (COL_JET): - { - hsl_to_rgb_jet(h, s, l, rgb); - break; - } - case (COL_HSLUV): - { - hsluv2rgb(h, 100.0*s, l, &r, &g, &b); - rgb[0] = r*100.0; - rgb[1] = g*100.0; - rgb[2] = b*100.0; - break; - } - } -// if (HSLUV_COLORS) -// { -// hsluv2rgb(h, 100.0*s, l, &r, &g, &b); -// rgb[0] = r*100.0; -// rgb[1] = g*100.0; -// rgb[2] = b*100.0; -// } -// else hsl_to_rgb_jet(h, s, l, rgb); -} - -double color_amplitude(double value, double scale, int time) -/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */ -{ - return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION))); -} -void color_scheme(int scheme, double value, double scale, int time, double rgb[3]) /* color scheme */ -{ - double hue, y, r, amplitude; - int intpart; - - /* saturation = r, luminosity = y */ - switch (scheme) { - case (C_LUM): - { - hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; - if (hue < 0.0) hue += 360.0; - if (hue >= 360.0) hue -= 360.0; - r = 0.9; - amplitude = color_amplitude(value, scale, time); - y = LUMMEAN + amplitude*LUMAMP; - intpart = (int)y; - y -= (double)intpart; - hsl_to_rgb(hue, r, y, rgb); - break; - } - case (C_HUE): - { - r = 0.9; - amplitude = color_amplitude(value, scale, time); - y = 0.5; - hue = HUEMEAN + amplitude*HUEAMP; - if (hue < 0.0) hue += 360.0; - if (hue >= 360.0) hue -= 360.0; - hsl_to_rgb(hue, r, y, rgb); - break; - } - } -} - - -void color_scheme_lum(int scheme, double value, double scale, int time, double lum, double rgb[3]) /* color scheme */ -{ - double hue, y, r, amplitude; - int intpart; - - /* saturation = r, luminosity = y */ - switch (scheme) { - case (C_LUM): - { - hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; - if (hue < 0.0) hue += 360.0; - if (hue >= 360.0) hue -= 360.0; - r = 0.9; - amplitude = color_amplitude(value, scale, time); - y = LUMMEAN + amplitude*LUMAMP; - intpart = (int)y; - y -= (double)intpart; - hsl_to_rgb(hue, r, y, rgb); - break; - } - case (C_HUE): - { - r = 0.9; - amplitude = color_amplitude(value, scale, time); - y = lum; - hue = HUEMEAN + amplitude*HUEAMP; - if (hue < 0.0) hue += 360.0; - if (hue >= 360.0) hue -= 360.0; - hsl_to_rgb(hue, r, y, rgb); - break; - } - } -} void blank() @@ -502,7 +364,7 @@ void init_circle_config() /* for billiard shape D_CIRCLES */ { int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted; - double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU; + double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU, xx[4], yy[4]; short int active_poisson[NMAXCIRCLES], far; switch (CIRCLE_PATTERN) { @@ -621,6 +483,42 @@ void init_circle_config() } break; } + case (C_LASER): + { + ncircles = 17; + + xx[0] = 0.5*(X_SHOOTER + X_TARGET); + xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER); + xx[2] = -xx[0]; + xx[3] = -xx[1]; + + yy[0] = 0.5*(Y_SHOOTER + Y_TARGET); + yy[1] = 1.0 - 0.5*(Y_TARGET - Y_SHOOTER); + yy[2] = -yy[0]; + yy[3] = -yy[1]; + + for (i = 0; i < 4; i++) + for (j = 0; j < 4; j++) + { + circlex[4*i + j] = xx[i]; + circley[4*i + j] = yy[j]; + + } + + circlex[ncircles - 1] = X_TARGET; + circley[ncircles - 1] = Y_TARGET; + + for (i=0; i= 0.0)&&(v2 >= 0.0)&&(v3 >= 0.0)) return(1); + else return(0); +} + + int xy_in_billiard(double x, double y) /* returns 1 if (x,y) represents a point in the billiard */ // double x, y; { double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width; int i, j, k, k1, k2, condition; + static double vertex[9][2], wertex[9][2]; + static int first = 1; switch (B_DOMAIN) { case (D_RECTANGLE): @@ -1048,6 +1039,33 @@ int xy_in_billiard(double x, double y) else return(1); } } + case (D_ISOSPECTRAL): + { + if (first) + { + compute_isospectral_coordinates(0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, vertex); + compute_isospectral_coordinates(1, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, wertex); + first = 0; + } + /* 1st triangle */ + condition = xy_in_triangle(x, y, vertex[0], vertex[1], vertex[2]); + condition += xy_in_triangle(x, y, vertex[0], vertex[4], vertex[1]); + condition += xy_in_triangle(x, y, vertex[1], vertex[5], vertex[2]); + condition += xy_in_triangle(x, y, vertex[0], vertex[2], vertex[3]); + condition += xy_in_triangle(x, y, vertex[1], vertex[4], vertex[7]); + condition += xy_in_triangle(x, y, vertex[2], vertex[5], vertex[8]); + condition += xy_in_triangle(x, y, vertex[0], vertex[3], vertex[6]); + + /* 2nd triangle */ + condition += xy_in_triangle(x, y, wertex[0], wertex[1], wertex[2]); + condition += xy_in_triangle(x, y, wertex[0], wertex[4], wertex[1]); + condition += xy_in_triangle(x, y, wertex[1], wertex[5], wertex[2]); + condition += xy_in_triangle(x, y, wertex[0], wertex[2], wertex[3]); + condition += xy_in_triangle(x, y, wertex[0], wertex[7], wertex[4]); + condition += xy_in_triangle(x, y, wertex[1], wertex[8], wertex[5]); + condition += xy_in_triangle(x, y, wertex[2], wertex[6], wertex[3]); + return(condition >= 1); + } case (D_CIRCLES): { for (i = 0; i < ncircles; i++) @@ -1060,6 +1078,19 @@ int xy_in_billiard(double x, double y) } return(1); } + case (D_CIRCLES_IN_RECT): /* returns 2 inside circles, 0 outside rectangle */ + { + for (i = 0; i < ncircles; i++) + if (circleactive[i]) + { + x1 = circlex[i]; + y1 = circley[i]; + r2 = circlerad[i]*circlerad[i]; + if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(2); + } + if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0); + else return(1); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -1225,10 +1256,22 @@ void toka_lineto(double x1, double y1) glVertex2d(pos[0], pos[1]); } +void iso_lineto(double z[2]) +/* draws boundary segments of isospectral billiard */ +{ + double pos[2]; + + xy_to_pos(z[0], z[1], pos); + glVertex2d(pos[0], pos[1]); +} + + void draw_billiard() /* draws the billiard boundary */ { double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax; + static double vertex[9][2], wertex[9][2]; int i, j, k, k1, k2, mr2; + static int first = 1; if (BLACK) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); @@ -1635,7 +1678,7 @@ void draw_billiard() /* draws the billiard boundary */ glColor3f(0.3, 0.3, 0.3); for (k=0; k= INITIAL_TIME)&&(DOUBLE_MOVIE)) { - draw_wave(phi, psi, xy_in, scale, i, PLOT_B); +// draw_wave(phi, psi, xy_in, scale, i, PLOT_B); + draw_wave_e(phi, psi, total_energy, xy_in, scale, i, PLOT_B); + if (DRAW_COLOR_SCHEME) draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, PLOT_B, -12.0, 12.0); draw_billiard(); glutSwapBuffers(); save_frame_counter(NSTEPS + 21 + counter); @@ -434,14 +491,16 @@ void animation() { if (DOUBLE_MOVIE) { - draw_wave(phi, psi, xy_in, scale, i, PLOT); +// draw_wave(phi, psi, xy_in, scale, i, PLOT); + draw_wave_e(phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT); draw_billiard(); glutSwapBuffers(); } for (i=0; i 0 and x < 0 do not communicate - phi is wave height, psi is phi at time t-1 */ +{ + int i, j; + double xy[2], dist2; + + for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); + else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); break; } case (P_MIXED): @@ -352,7 +377,9 @@ void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], sho { energy = compute_energy(phi, psi, xy_in, i, j); total_energy[i][j] += energy; - color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme(COLOR_SCHEME, total_energy[i][j]/(double)(time+1), scale, time, rgb); break; } } @@ -369,6 +396,59 @@ void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], sho } +void draw_color_scheme(double x1, double y1, double x2, double y2, int plot, double min, double max) +{ + int j, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; + double y, dy, dy_e, rgb[3], value; + + xy_to_ij(x1, y1, ij_botleft); + xy_to_ij(x2, y2, ij_topright); + + imin = ij_botleft[0]; + imax = ij_topright[0]; + jmin = ij_botleft[1]; + jmax = ij_topright[1]; + + glBegin(GL_QUADS); + dy = (max - min)/((double)(jmax - jmin)); + dy_e = max/((double)(jmax - jmin)); + + for (j = jmin; j < jmax; j++) + { + switch (plot) { + case (P_AMPLITUDE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + case (P_ENERGY): + { + value = dy_e*(double)(j - jmin); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); + else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + case (P_MEAN_ENERGY): + { + value = dy_e*(double)(j - jmin); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); + else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + } + glColor3f(rgb[0], rgb[1], rgb[2]); + glVertex2i(imin, j); + glVertex2i(imax, j); + glVertex2i(imax, j+1); + glVertex2i(imin, j+1); + } + glEnd (); + + glColor3f(1.0, 1.0, 1.0); + draw_rectangle(x1, y1, x2, y2); +} + diff --git a/wave_comparison.c b/wave_comparison.c index 5df0673..6539b91 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -81,6 +81,17 @@ #define NGRIDX 20 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ @@ -133,9 +144,11 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.2 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.002 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ +#define INITIAL_AMP 0.75 /* amplitude of initial condition */ +// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ +#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ @@ -143,7 +156,7 @@ /* Color schemes */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -166,8 +179,6 @@ #define VMAX 5.0 /* max value of wave amplitude */ -#include "hsluv.c" - #include "global_pdes.c" /* constants and global variables */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ diff --git a/wave_energy.c b/wave_energy.c index 572d0f6..af88bac 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -50,23 +50,23 @@ #define NX 1280 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ -#define XMIN -1.777777778 -#define XMAX 1.777777778 /* x interval */ -#define YMIN -1.0 -#define YMAX 1.0 /* y interval for 9/16 aspect ratio */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +// #define XMIN -1.777777778 +// #define XMAX 1.777777778 /* x interval */ +// #define YMIN -1.0 +// #define YMAX 1.0 /* y interval for 9/16 aspect ratio */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 15 /* choice of domain shape, see list in global_pdes.c */ -#define B_DOMAIN_B 15 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 2 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 11 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ @@ -85,6 +85,17 @@ #define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +#define ISO_XSHIFT_LEFT -1.65 +#define ISO_XSHIFT_RIGHT 0.4 +#define ISO_YSHIFT_LEFT -0.05 +#define ISO_YSHIFT_RIGHT -0.05 +#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ + /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ @@ -98,14 +109,8 @@ #define AMPLITUDE 0.025 /* amplitude of periodic excitation */ #define COURANT 0.02 /* Courant number */ #define COURANTB 0.004 /* Courant number in medium B */ -// #define COURANTB 0.005 /* Courant number in medium B */ -// #define COURANTB 0.008 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -// #define GAMMA 1.0e-8 /* damping factor in wave equation */ #define GAMMAB 1.0e-8 /* damping factor in wave equation */ -// #define GAMMAB 1.0e-6 /* damping factor in wave equation */ -// #define GAMMAB 2.0e-4 /* damping factor in wave equation */ -// #define GAMMAB 2.5e-4 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -122,7 +127,7 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 4500 /* number of frames of movie */ +#define NSTEPS 3750 /* number of frames of movie */ #define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 200 /* time after which to start saving frames */ @@ -170,8 +175,6 @@ #define VMAX 5.0 /* max value of wave amplitude */ -#include "hsluv.c" - #include "global_pdes.c" /* constants and global variables */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */