A Comparative Study of Limiting Strategies in Discontinuous Galerkin Schemes for the $M_1$ Model of Radiation Transport
The $M_1$ minimum entropy moment system is a system of hyperbolic balance laws that approximates the radiation transport equation, and has many desirable properties. Among them are symmetric hyperbolicity, entropy decay, moment realizability, and correct behavior in the diffusion and free-streaming limits. However, numerical difficulties arise when approximating the solution of the $M_1$ model by high order numerical schemes; namely maintaining the realizability of the numerical solution and controlling spurious oscillations. In this paper, we extend a previously constructed one-dimensional realizability limiting strategy to 2D. In addition, we perform a numerical study of various combinations of the realizability limiter and the TVBM local slope limiter on a third order Discontinuous Galerkin (DG) scheme on both triangular and rectangular meshes. In several test cases, we demonstrate that in general, a combination of the realizability limiter and a TVBM limiter is necessary to obtain a robust and accurate numerical scheme. Our code is published so that all results can be reproduced by the reader.
colormap.tex
\def\mycolormap{                          
rgb(0pt)=(0.2081,0.1663,0.5292);          
rgb(1pt)=(0.208355,0.16778,0.532238);     
rgb(2pt)=(0.208611,0.169261,0.535275);    
rgb(3pt)=(0.208866,0.170741,0.538313);    
rgb(4pt)=(0.209121,0.172222,0.54135);     
rgb(5pt)=(0.209376,0.173702,0.544388);    
rgb(6pt)=(0.209632,0.175183,0.547425);    
rgb(7pt)=(0.209887,0.176663,0.550463);    
rgb(8pt)=(0.210134,0.178144,0.553505);    
rgb(9pt)=(0.210338,0.179624,0.556568);    
rgb(10pt)=(0.210542,0.181105,0.559631);   
rgb(11pt)=(0.210746,0.182585,0.562694);   
rgb(12pt)=(0.210944,0.184066,0.565763);   
rgb(13pt)=(0.211123,0.185546,0.568852);   
rgb(14pt)=(0.211302,0.187027,0.57194);    
rgb(15pt)=(0.21148,0.188507,0.575029);    
rgb(16pt)=(0.211642,0.189996,0.578117);   
rgb(17pt)=(0.21177,0.191502,0.581206);    
rgb(18pt)=(0.211897,0.193008,0.584295);   
rgb(19pt)=(0.212025,0.194514,0.587383);   
rgb(20pt)=(0.212132,0.19602,0.590472);    
rgb(21pt)=(0.212208,0.197526,0.59356);    
rgb(22pt)=(0.212285,0.199032,0.596649);   
rgb(23pt)=(0.212361,0.200538,0.599738);   
rgb(24pt)=(0.212413,0.202044,0.602839);   
rgb(25pt)=(0.212438,0.20355,0.605953);    
rgb(26pt)=(0.212464,0.205056,0.609067);   
rgb(27pt)=(0.212489,0.206562,0.612181);   
rgb(28pt)=(0.212471,0.208083,0.61531);    
rgb(29pt)=(0.21242,0.209614,0.61845);     
rgb(30pt)=(0.212368,0.211146,0.621589);   
rgb(31pt)=(0.212317,0.212677,0.624729);   
rgb(32pt)=(0.212216,0.214209,0.627868);   
rgb(33pt)=(0.212088,0.215741,0.631008);   
rgb(34pt)=(0.211961,0.217272,0.634148);   
rgb(35pt)=(0.211833,0.218804,0.637287);   
rgb(36pt)=(0.211668,0.220354,0.640446);   
rgb(37pt)=(0.211489,0.221911,0.643611);   
rgb(38pt)=(0.21131,0.223468,0.646776);    
rgb(39pt)=(0.211132,0.225025,0.649941);   
rgb(40pt)=(0.210848,0.226603,0.653107);   
rgb(41pt)=(0.210541,0.228186,0.656272);   
rgb(42pt)=(0.210235,0.229768,0.659437);   
rgb(43pt)=(0.209929,0.231351,0.662602);   
rgb(44pt)=(0.209553,0.232934,0.665767);   
rgb(45pt)=(0.20917,0.234516,0.668932);    
rgb(46pt)=(0.208787,0.236099,0.672098);   
rgb(47pt)=(0.208405,0.237681,0.675263);   
rgb(48pt)=(0.20787,0.239289,0.678453);    
rgb(49pt)=(0.207334,0.240897,0.681644);   
rgb(50pt)=(0.206798,0.242505,0.684835);   
rgb(51pt)=(0.206255,0.244114,0.688025);   
rgb(52pt)=(0.205617,0.245722,0.691216);   
rgb(53pt)=(0.204979,0.24733,0.694407);    
rgb(54pt)=(0.204341,0.248938,0.697597);   
rgb(55pt)=(0.203675,0.250554,0.700792);   
rgb(56pt)=(0.202858,0.252213,0.704008);   
rgb(57pt)=(0.202041,0.253872,0.707224);   
rgb(58pt)=(0.201225,0.255531,0.710441);   
rgb(59pt)=(0.200372,0.257184,0.713657);   
rgb(60pt)=(0.199402,0.258818,0.716873);   
rgb(61pt)=(0.198432,0.260452,0.720089);   
rgb(62pt)=(0.197462,0.262085,0.723305);   
rgb(63pt)=(0.196419,0.263735,0.726522);   
rgb(64pt)=(0.195219,0.26542,0.729738);    
rgb(65pt)=(0.19402,0.267105,0.732954);    
rgb(66pt)=(0.19282,0.268789,0.73617);     
rgb(67pt)=(0.191549,0.270474,0.739386);   
rgb(68pt)=(0.19017,0.272159,0.742603);    
rgb(69pt)=(0.188792,0.273843,0.745819);   
rgb(70pt)=(0.187414,0.275528,0.749035);   
rgb(71pt)=(0.1859,0.277237,0.752264);     
rgb(72pt)=(0.184241,0.278973,0.755505);   
rgb(73pt)=(0.182581,0.280709,0.758747);   
rgb(74pt)=(0.180922,0.282444,0.761989);   
rgb(75pt)=(0.179133,0.284209,0.765245);   
rgb(76pt)=(0.177244,0.285996,0.768512);   
rgb(77pt)=(0.175356,0.287783,0.77178);    
rgb(78pt)=(0.173467,0.289569,0.775047);   
rgb(79pt)=(0.171363,0.291406,0.778314);   
rgb(80pt)=(0.169142,0.293269,0.781581);   
rgb(81pt)=(0.166922,0.295132,0.784849);   
rgb(82pt)=(0.164701,0.296996,0.788116);   
rgb(83pt)=(0.162238,0.298934,0.791365);   
rgb(84pt)=(0.159686,0.300899,0.794606);   
rgb(85pt)=(0.157133,0.302865,0.797848);   
rgb(86pt)=(0.15458,0.30483,0.80109);      
rgb(87pt)=(0.151738,0.306858,0.804352);   
rgb(88pt)=(0.148828,0.3089,0.80762);      
rgb(89pt)=(0.145918,0.310942,0.810887);   
rgb(90pt)=(0.143008,0.312984,0.814154);   
rgb(91pt)=(0.139687,0.31514,0.81733);     
rgb(92pt)=(0.136318,0.31731,0.820495);    
rgb(93pt)=(0.132949,0.319479,0.82366);    
rgb(94pt)=(0.129579,0.321649,0.826826);   
rgb(95pt)=(0.125811,0.323918,0.829841);   
rgb(96pt)=(0.122033,0.32619,0.832853);    
rgb(97pt)=(0.118256,0.328462,0.835865);   
rgb(98pt)=(0.114458,0.330737,0.838862);   
rgb(99pt)=(0.110349,0.333059,0.841619);   
rgb(100pt)=(0.106239,0.335382,0.844376);  
rgb(101pt)=(0.102129,0.337705,0.847132);  
rgb(102pt)=(0.0979874,0.340021,0.849835); 
rgb(103pt)=(0.093648,0.342292,0.852209);  
rgb(104pt)=(0.0893087,0.344564,0.854583); 
rgb(105pt)=(0.0849694,0.346836,0.856957); 
rgb(106pt)=(0.08063,0.349091,0.859234);   
rgb(107pt)=(0.0762907,0.351286,0.861174); 
rgb(108pt)=(0.0719514,0.353481,0.863114); 
rgb(109pt)=(0.067612,0.355676,0.865053);  
rgb(110pt)=(0.0633195,0.357817,0.866853); 
rgb(111pt)=(0.0591333,0.359833,0.868333); 
rgb(112pt)=(0.0549471,0.36185,0.869814);  
rgb(113pt)=(0.050761,0.363866,0.871294);  
rgb(114pt)=(0.0466838,0.365823,0.872626); 
rgb(115pt)=(0.0427784,0.367687,0.873724); 
rgb(116pt)=(0.038873,0.36955,0.874821);   
rgb(117pt)=(0.0349676,0.371414,0.875919); 
rgb(118pt)=(0.0315066,0.373217,0.876872); 
rgb(119pt)=(0.0285456,0.374953,0.877664); 
rgb(120pt)=(0.0255847,0.376688,0.878455); 
rgb(121pt)=(0.0226237,0.378424,0.879246); 
rgb(122pt)=(0.0202132,0.380061,0.879868); 
rgb(123pt)=(0.0182477,0.381618,0.880353); 
rgb(124pt)=(0.0162823,0.383175,0.880838); 
rgb(125pt)=(0.0143168,0.384732,0.881323); 
rgb(126pt)=(0.0127892,0.386241,0.881695); 
rgb(127pt)=(0.0115129,0.387721,0.882001); 
rgb(128pt)=(0.0102366,0.389202,0.882307); 
rgb(129pt)=(0.00896036,0.390682,0.882614);
rgb(130pt)=(0.00812372,0.392089,0.88281); 
rgb(131pt)=(0.00746006,0.393468,0.882963);
rgb(132pt)=(0.0067964,0.394846,0.883116); 
rgb(133pt)=(0.00613273,0.396224,0.883269);
rgb(134pt)=(0.00581622,0.397562,0.88332); 
rgb(135pt)=(0.00558649,0.398889,0.883346);
rgb(136pt)=(0.00535676,0.400217,0.883371);
rgb(137pt)=(0.00512703,0.401544,0.883397);
rgb(138pt)=(0.00516757,0.402804,0.883332);
rgb(139pt)=(0.00524414,0.404054,0.883256);
rgb(140pt)=(0.00532072,0.405305,0.883179);
rgb(141pt)=(0.0053973,0.406556,0.883103); 
rgb(142pt)=(0.00572012,0.407757,0.882952);
rgb(143pt)=(0.00605195,0.408957,0.882799);
rgb(144pt)=(0.00638378,0.410157,0.882646);
rgb(145pt)=(0.00672643,0.411355,0.882489);
rgb(146pt)=(0.00728799,0.412529,0.882259);
rgb(147pt)=(0.00784955,0.413704,0.88203); 
rgb(148pt)=(0.00841111,0.414878,0.8818);  
rgb(149pt)=(0.00898919,0.416045,0.881564);
rgb(150pt)=(0.00967838,0.417168,0.881283);
rgb(151pt)=(0.0103676,0.418292,0.881002); 
rgb(152pt)=(0.0110568,0.419415,0.880721); 
rgb(153pt)=(0.011773,0.420532,0.880435);  
rgb(154pt)=(0.0125898,0.42163,0.880129);  
rgb(155pt)=(0.0134066,0.422728,0.879823); 
rgb(156pt)=(0.0142234,0.423825,0.879516); 
rgb(157pt)=(0.0150703,0.424915,0.879195); 
rgb(158pt)=(0.0159892,0.425987,0.878838); 
rgb(159pt)=(0.0169081,0.427059,0.87848);  
rgb(160pt)=(0.017827,0.428132,0.878123);  
rgb(161pt)=(0.0187748,0.429194,0.877746); 
rgb(162pt)=(0.0197703,0.430241,0.877338); 
rgb(163pt)=(0.0207658,0.431287,0.876929); 
rgb(164pt)=(0.0217613,0.432334,0.876521); 
rgb(165pt)=(0.0227802,0.43338,0.876113);  
rgb(166pt)=(0.0238267,0.434427,0.875704); 
rgb(167pt)=(0.0248733,0.435473,0.875296); 
rgb(168pt)=(0.0259198,0.43652,0.874887);  
rgb(169pt)=(0.0269802,0.437553,0.874451); 
rgb(170pt)=(0.0280523,0.438574,0.873992); 
rgb(171pt)=(0.0291243,0.439595,0.873532); 
rgb(172pt)=(0.0301964,0.440616,0.873073); 
rgb(173pt)=(0.0312844,0.441621,0.872614); 
rgb(174pt)=(0.032382,0.442616,0.872154);  
rgb(175pt)=(0.0334796,0.443612,0.871695); 
rgb(176pt)=(0.0345772,0.444607,0.871235); 
rgb(177pt)=(0.0357108,0.445603,0.870758); 
rgb(178pt)=(0.0368595,0.446598,0.870273); 
rgb(179pt)=(0.0380081,0.447594,0.869788); 
rgb(180pt)=(0.0391568,0.448589,0.869303); 
rgb(181pt)=(0.0402652,0.449565,0.868798); 
rgb(182pt)=(0.0413628,0.450535,0.868287); 
rgb(183pt)=(0.0424604,0.451505,0.867777); 
rgb(184pt)=(0.043558,0.452474,0.867266);  
rgb(185pt)=(0.0445889,0.453444,0.866756); 
rgb(186pt)=(0.0456099,0.454414,0.866245); 
rgb(187pt)=(0.0466309,0.455384,0.865735); 
rgb(188pt)=(0.047652,0.456354,0.865224);  
rgb(189pt)=(0.0486,0.457324,0.864714);    
rgb(190pt)=(0.0495444,0.458294,0.864203); 
rgb(191pt)=(0.0504889,0.459264,0.863692); 
rgb(192pt)=(0.0514315,0.460234,0.863181); 
rgb(193pt)=(0.0523249,0.461204,0.862645); 
rgb(194pt)=(0.0532183,0.462174,0.862109); 
rgb(195pt)=(0.0541117,0.463144,0.861573); 
rgb(196pt)=(0.0549991,0.464111,0.861034); 
rgb(197pt)=(0.0558414,0.465056,0.860472); 
rgb(198pt)=(0.0566838,0.466,0.859911);    
rgb(199pt)=(0.0575261,0.466944,0.859349); 
rgb(200pt)=(0.0583532,0.467889,0.858793); 
rgb(201pt)=(0.0591189,0.468833,0.858257); 
rgb(202pt)=(0.0598847,0.469778,0.857721); 
rgb(203pt)=(0.0606505,0.470722,0.857185); 
rgb(204pt)=(0.0614018,0.471667,0.856641); 
rgb(205pt)=(0.0621165,0.472611,0.85608);  
rgb(206pt)=(0.0628312,0.473556,0.855518); 
rgb(207pt)=(0.0635459,0.4745,0.854957);   
rgb(208pt)=(0.064242,0.475444,0.854405);  
rgb(209pt)=(0.0649057,0.476389,0.853868); 
rgb(210pt)=(0.0655694,0.477333,0.853332); 
rgb(211pt)=(0.066233,0.478278,0.852796);  
rgb(212pt)=(0.0668625,0.479222,0.852249); 
rgb(213pt)=(0.0674495,0.480167,0.851687); 
rgb(214pt)=(0.0680366,0.481111,0.851126); 
rgb(215pt)=(0.0686237,0.482056,0.850564); 
rgb(216pt)=(0.0691838,0.483,0.850003);    
rgb(217pt)=(0.0697198,0.483944,0.849441); 
rgb(218pt)=(0.0702559,0.484889,0.84888);  
rgb(219pt)=(0.0707919,0.485833,0.848318); 
rgb(220pt)=(0.0712967,0.486778,0.847772); 
rgb(221pt)=(0.0717817,0.487722,0.847236); 
rgb(222pt)=(0.0722667,0.488667,0.8467);   
rgb(223pt)=(0.0727517,0.489611,0.846164); 
rgb(224pt)=(0.0732012,0.490573,0.845628); 
rgb(225pt)=(0.0736351,0.491543,0.845092); 
rgb(226pt)=(0.0740691,0.492513,0.844556); 
rgb(227pt)=(0.074503,0.493483,0.84402);   
rgb(228pt)=(0.0748973,0.494433,0.843484); 
rgb(229pt)=(0.0752802,0.495378,0.842948); 
rgb(230pt)=(0.0756631,0.496322,0.842412); 
rgb(231pt)=(0.0760459,0.497267,0.841876); 
rgb(232pt)=(0.0763631,0.498233,0.841362); 
rgb(233pt)=(0.0766694,0.499203,0.840851); 
rgb(234pt)=(0.0769757,0.500173,0.840341); 
rgb(235pt)=(0.077282,0.501143,0.83983);   
rgb(236pt)=(0.0775162,0.502137,0.83932);  
rgb(237pt)=(0.0777459,0.503132,0.838809); 
rgb(238pt)=(0.0779757,0.504128,0.838298); 
rgb(239pt)=(0.0782042,0.505123,0.837789); 
rgb(240pt)=(0.0783829,0.506093,0.837304); 
rgb(241pt)=(0.0785616,0.507063,0.836819); 
rgb(242pt)=(0.0787402,0.508033,0.836334); 
rgb(243pt)=(0.0789135,0.509008,0.835851); 
rgb(244pt)=(0.0790411,0.510029,0.835392); 
rgb(245pt)=(0.0791688,0.51105,0.834932);  
rgb(246pt)=(0.0792964,0.512071,0.834473); 
rgb(247pt)=(0.0794048,0.513092,0.834018); 
rgb(248pt)=(0.0794303,0.514113,0.833584); 
rgb(249pt)=(0.0794559,0.515134,0.83315);  
rgb(250pt)=(0.0794814,0.516155,0.832717); 
rgb(251pt)=(0.0794862,0.517183,0.832289); 
rgb(252pt)=(0.0794351,0.51823,0.831881);  
rgb(253pt)=(0.0793841,0.519276,0.831473); 
rgb(254pt)=(0.079333,0.520323,0.831064);  
rgb(255pt)=(0.079255,0.521369,0.830665);  
rgb(256pt)=(0.0791273,0.522416,0.830282); 
rgb(257pt)=(0.0789997,0.523462,0.829899); 
rgb(258pt)=(0.0788721,0.524509,0.829516); 
rgb(259pt)=(0.0786889,0.525589,0.829156); 
rgb(260pt)=(0.0784336,0.526712,0.828824); 
rgb(261pt)=(0.0781784,0.527835,0.828492); 
rgb(262pt)=(0.0779231,0.528958,0.82816);  
rgb(263pt)=(0.077615,0.530081,0.827868);  
rgb(264pt)=(0.0772577,0.531205,0.827613); 
rgb(265pt)=(0.0769003,0.532328,0.827357); 
rgb(266pt)=(0.0765429,0.533451,0.827102); 
rgb(267pt)=(0.0761243,0.534589,0.826862); 
rgb(268pt)=(0.0756649,0.535738,0.826632); 
rgb(269pt)=(0.0752054,0.536886,0.826403); 
rgb(270pt)=(0.0747459,0.538035,0.826173); 
rgb(271pt)=(0.0742168,0.539219,0.825961); 
rgb(272pt)=(0.0736553,0.540418,0.825756); 
rgb(273pt)=(0.0730937,0.541618,0.825552); 
rgb(274pt)=(0.0725321,0.542818,0.825348); 
rgb(275pt)=(0.0718925,0.544037,0.825183); 
rgb(276pt)=(0.0712288,0.545262,0.82503);  
rgb(277pt)=(0.0705652,0.546487,0.824877); 
rgb(278pt)=(0.0699015,0.547713,0.824723); 
rgb(279pt)=(0.0691514,0.548938,0.824614); 
rgb(280pt)=(0.0683856,0.550163,0.824511); 
rgb(281pt)=(0.0676198,0.551388,0.824409); 
rgb(282pt)=(0.0668541,0.552614,0.824307); 
rgb(283pt)=(0.0660408,0.553886,0.824205); 
rgb(284pt)=(0.065224,0.555162,0.824103);  
rgb(285pt)=(0.0644072,0.556439,0.824001); 
rgb(286pt)=(0.0635892,0.557715,0.823899); 
rgb(287pt)=(0.0626703,0.558991,0.823848); 
rgb(288pt)=(0.0617514,0.560268,0.823797); 
rgb(289pt)=(0.0608324,0.561544,0.823746); 
rgb(290pt)=(0.0599087,0.56282,0.823693);  
rgb(291pt)=(0.0589387,0.564096,0.823616); 
rgb(292pt)=(0.0579688,0.565373,0.82354);  
rgb(293pt)=(0.0569988,0.566649,0.823463); 
rgb(294pt)=(0.0560243,0.567925,0.823386); 
rgb(295pt)=(0.0550288,0.569202,0.82331);  
rgb(296pt)=(0.0540333,0.570478,0.823233); 
rgb(297pt)=(0.0530378,0.571754,0.823157); 
rgb(298pt)=(0.0520423,0.57303,0.82308);   
rgb(299pt)=(0.0510468,0.574307,0.823004); 
rgb(300pt)=(0.0500514,0.575583,0.822927); 
rgb(301pt)=(0.0490559,0.576859,0.82285);  
rgb(302pt)=(0.0480604,0.578127,0.822756); 
rgb(303pt)=(0.0470649,0.579377,0.822629); 
rgb(304pt)=(0.0460694,0.580628,0.822501); 
rgb(305pt)=(0.0450739,0.581879,0.822374); 
rgb(306pt)=(0.0441,0.583119,0.822235);    
rgb(307pt)=(0.0431556,0.584344,0.822082); 
rgb(308pt)=(0.0422111,0.585569,0.821929); 
rgb(309pt)=(0.0412667,0.586795,0.821776); 
rgb(310pt)=(0.0403351,0.58802,0.821597);  
rgb(311pt)=(0.0394162,0.589245,0.821392); 
rgb(312pt)=(0.0384973,0.59047,0.821188);  
rgb(313pt)=(0.0375784,0.591695,0.820984); 
rgb(314pt)=(0.0367495,0.592891,0.820735); 
rgb(315pt)=(0.0359838,0.594065,0.820454); 
rgb(316pt)=(0.035218,0.595239,0.820173);  
rgb(317pt)=(0.0344523,0.596413,0.819892); 
rgb(318pt)=(0.0337721,0.597553,0.819595); 
rgb(319pt)=(0.0331339,0.598676,0.819288); 
rgb(320pt)=(0.0324958,0.599799,0.818982); 
rgb(321pt)=(0.0318577,0.600923,0.818676); 
rgb(322pt)=(0.0312964,0.602026,0.818312); 
rgb(323pt)=(0.0307604,0.603124,0.817929); 
rgb(324pt)=(0.0302243,0.604222,0.817546); 
rgb(325pt)=(0.0296883,0.605319,0.817163); 
rgb(326pt)=(0.0292375,0.606395,0.816738); 
rgb(327pt)=(0.0288036,0.607468,0.816304); 
rgb(328pt)=(0.0283697,0.60854,0.81587);   
rgb(329pt)=(0.0279357,0.609612,0.815436); 
rgb(330pt)=(0.0275721,0.610637,0.814955); 
rgb(331pt)=(0.0272147,0.611658,0.81447);  
rgb(332pt)=(0.0268574,0.612679,0.813985); 
rgb(333pt)=(0.0265,0.6137,0.8135);        
rgb(334pt)=(0.0262447,0.614695,0.812964); 
rgb(335pt)=(0.0259895,0.615691,0.812428); 
rgb(336pt)=(0.0257342,0.616686,0.811892); 
rgb(337pt)=(0.0254853,0.61768,0.811352);  
rgb(338pt)=(0.0253066,0.61865,0.810765);  
rgb(339pt)=(0.0251279,0.61962,0.810177);  
rgb(340pt)=(0.0249492,0.62059,0.80959);   
rgb(341pt)=(0.024779,0.621551,0.808995);  
rgb(342pt)=(0.0246514,0.62247,0.808357);  
rgb(343pt)=(0.0245237,0.623389,0.807719); 
rgb(344pt)=(0.0243961,0.624308,0.80708);  
rgb(345pt)=(0.0242748,0.625221,0.80643);  
rgb(346pt)=(0.0241727,0.626114,0.805741); 
rgb(347pt)=(0.0240706,0.627008,0.805051); 
rgb(348pt)=(0.0239685,0.627901,0.804362); 
rgb(349pt)=(0.0238832,0.628786,0.803656); 
rgb(350pt)=(0.0238321,0.629654,0.802916); 
rgb(351pt)=(0.0237811,0.630522,0.802176); 
rgb(352pt)=(0.02373,0.631389,0.801435);   
rgb(353pt)=(0.023679,0.632247,0.800685);  
rgb(354pt)=(0.0236279,0.633089,0.799919); 
rgb(355pt)=(0.0235769,0.633932,0.799153); 
rgb(356pt)=(0.0235258,0.634774,0.798387); 
rgb(357pt)=(0.0234748,0.635604,0.797596); 
rgb(358pt)=(0.0234237,0.63642,0.79678);   
rgb(359pt)=(0.0233727,0.637237,0.795963); 
rgb(360pt)=(0.0233216,0.638054,0.795146); 
rgb(361pt)=(0.0232706,0.638856,0.794329); 
rgb(362pt)=(0.0232195,0.639647,0.793512); 
rgb(363pt)=(0.0231685,0.640439,0.792695); 
rgb(364pt)=(0.0231174,0.64123,0.791879);  
rgb(365pt)=(0.0230832,0.642005,0.791011); 
rgb(366pt)=(0.0230577,0.64277,0.790118);  
rgb(367pt)=(0.0230321,0.643536,0.789225); 
rgb(368pt)=(0.0230066,0.644302,0.788331); 
rgb(369pt)=(0.0229811,0.645049,0.787438); 
rgb(370pt)=(0.0229556,0.645789,0.786544); 
rgb(371pt)=(0.02293,0.646529,0.785651);   
rgb(372pt)=(0.0229045,0.647269,0.784758); 
rgb(373pt)=(0.022858,0.64801,0.783843);   
rgb(374pt)=(0.0228069,0.64875,0.782924);  
rgb(375pt)=(0.0227559,0.64949,0.782005);  
rgb(376pt)=(0.0227048,0.65023,0.781086);  
rgb(377pt)=(0.0227,0.650947,0.780144);    
rgb(378pt)=(0.0227,0.651662,0.7792);      
rgb(379pt)=(0.0227,0.652377,0.778256);    
rgb(380pt)=(0.0227,0.653092,0.777311);    
rgb(381pt)=(0.0228261,0.653781,0.776341); 
rgb(382pt)=(0.0229538,0.65447,0.775371);  
rgb(383pt)=(0.0230814,0.655159,0.774402); 
rgb(384pt)=(0.0232108,0.655849,0.77343);  
rgb(385pt)=(0.023364,0.656538,0.772434);  
rgb(386pt)=(0.0235171,0.657227,0.771439); 
rgb(387pt)=(0.0236703,0.657916,0.770443); 
rgb(388pt)=(0.0238312,0.658602,0.769444); 
rgb(389pt)=(0.0240354,0.659265,0.768423); 
rgb(390pt)=(0.0242396,0.659929,0.767402); 
rgb(391pt)=(0.0244438,0.660592,0.766381); 
rgb(392pt)=(0.0247021,0.661256,0.765354); 
rgb(393pt)=(0.025136,0.66192,0.764307);   
rgb(394pt)=(0.02557,0.662583,0.763261);   
rgb(395pt)=(0.0260039,0.663247,0.762214); 
rgb(396pt)=(0.0264541,0.663911,0.761168); 
rgb(397pt)=(0.026939,0.664574,0.760121);  
rgb(398pt)=(0.027424,0.665238,0.759074);  
rgb(399pt)=(0.027909,0.665902,0.758028);  
rgb(400pt)=(0.028445,0.666555,0.756971);  
rgb(401pt)=(0.0290577,0.667193,0.755899); 
rgb(402pt)=(0.0296703,0.667832,0.754827); 
rgb(403pt)=(0.0302829,0.66847,0.753755);  
rgb(404pt)=(0.030994,0.669095,0.752683);  
rgb(405pt)=(0.0318108,0.669708,0.751611); 
rgb(406pt)=(0.0326276,0.670321,0.750539); 
rgb(407pt)=(0.0334444,0.670933,0.749467); 
rgb(408pt)=(0.0343045,0.67156,0.748366);  
rgb(409pt)=(0.0351979,0.672198,0.747243); 
rgb(410pt)=(0.0360913,0.672837,0.74612);  
rgb(411pt)=(0.0369847,0.673475,0.744996); 
rgb(412pt)=(0.0380432,0.674096,0.743873); 
rgb(413pt)=(0.0391919,0.674709,0.74275);  
rgb(414pt)=(0.0403405,0.675322,0.741627); 
rgb(415pt)=(0.0414892,0.675934,0.740504); 
rgb(416pt)=(0.0427123,0.676528,0.739381); 
rgb(417pt)=(0.0439631,0.677115,0.738258); 
rgb(418pt)=(0.0452138,0.677702,0.737135); 
rgb(419pt)=(0.0464646,0.678289,0.736011); 
rgb(420pt)=(0.0477153,0.678897,0.734868); 
rgb(421pt)=(0.0489661,0.67951,0.733719);  
rgb(422pt)=(0.0502168,0.680123,0.73257);  
rgb(423pt)=(0.0514676,0.680735,0.731422); 
rgb(424pt)=(0.0529237,0.681325,0.73025);  
rgb(425pt)=(0.0544042,0.681912,0.729076); 
rgb(426pt)=(0.0558847,0.682499,0.727902); 
rgb(427pt)=(0.0573652,0.683086,0.726728); 
rgb(428pt)=(0.0587709,0.683673,0.725553); 
rgb(429pt)=(0.0601748,0.68426,0.724379);  
rgb(430pt)=(0.0615787,0.684847,0.723205); 
rgb(431pt)=(0.0629946,0.685435,0.722028); 
rgb(432pt)=(0.0646027,0.686022,0.720803); 
rgb(433pt)=(0.0662108,0.686609,0.719577); 
rgb(434pt)=(0.0678189,0.687196,0.718352); 
rgb(435pt)=(0.069427,0.687779,0.717131);  
rgb(436pt)=(0.0710351,0.688341,0.715931); 
rgb(437pt)=(0.0726432,0.688902,0.714731); 
rgb(438pt)=(0.0742514,0.689464,0.713532); 
rgb(439pt)=(0.0758709,0.690026,0.712326); 
rgb(440pt)=(0.07753,0.690587,0.711101);   
rgb(441pt)=(0.0791892,0.691149,0.709876); 
rgb(442pt)=(0.0808483,0.69171,0.70865);   
rgb(443pt)=(0.0825387,0.692272,0.707417); 
rgb(444pt)=(0.0843,0.692833,0.706167);    
rgb(445pt)=(0.0860613,0.693395,0.704916); 
rgb(446pt)=(0.0878225,0.693956,0.703665); 
rgb(447pt)=(0.089564,0.694518,0.702405);  
rgb(448pt)=(0.0912742,0.69508,0.701128);  
rgb(449pt)=(0.0929844,0.695641,0.699852); 
rgb(450pt)=(0.0946946,0.696203,0.698576); 
rgb(451pt)=(0.0965009,0.696752,0.697299); 
rgb(452pt)=(0.0984153,0.697288,0.696023); 
rgb(453pt)=(0.10033,0.697824,0.694747);   
rgb(454pt)=(0.102244,0.69836,0.693471);   
rgb(455pt)=(0.10413,0.698896,0.69218);    
rgb(456pt)=(0.105994,0.699432,0.690878);  
rgb(457pt)=(0.107857,0.699968,0.689577);  
rgb(458pt)=(0.10972,0.700505,0.688275);   
rgb(459pt)=(0.111632,0.701041,0.686973);  
rgb(460pt)=(0.113572,0.701577,0.685671);  
rgb(461pt)=(0.115512,0.702113,0.684369);  
rgb(462pt)=(0.117452,0.702649,0.683068);  
rgb(463pt)=(0.119429,0.703185,0.681747);  
rgb(464pt)=(0.12142,0.703721,0.68042);    
rgb(465pt)=(0.123411,0.704257,0.679093);  
rgb(466pt)=(0.125402,0.704793,0.677765);  
rgb(467pt)=(0.127372,0.705308,0.676438);  
rgb(468pt)=(0.129338,0.705819,0.675111);  
rgb(469pt)=(0.131303,0.706329,0.673783);  
rgb(470pt)=(0.133269,0.70684,0.672456);   
rgb(471pt)=(0.135369,0.70735,0.671084);   
rgb(472pt)=(0.137488,0.707861,0.669705);  
rgb(473pt)=(0.139607,0.708371,0.668327);  
rgb(474pt)=(0.141725,0.708882,0.666949);  
rgb(475pt)=(0.143795,0.709392,0.665595);  
rgb(476pt)=(0.145862,0.709903,0.664242);  
rgb(477pt)=(0.14793,0.710414,0.662889);   
rgb(478pt)=(0.150003,0.710924,0.661534);  
rgb(479pt)=(0.152198,0.711435,0.66013);   
rgb(480pt)=(0.154394,0.711945,0.658726);  
rgb(481pt)=(0.156589,0.712456,0.657322);  
rgb(482pt)=(0.158784,0.712963,0.655922);  
rgb(483pt)=(0.160979,0.713448,0.654543);  
rgb(484pt)=(0.163174,0.713933,0.653165);  
rgb(485pt)=(0.16537,0.714418,0.651786);   
rgb(486pt)=(0.16757,0.714908,0.650397);   
rgb(487pt)=(0.169791,0.715419,0.648968);  
rgb(488pt)=(0.172012,0.715929,0.647538);  
rgb(489pt)=(0.174232,0.71644,0.646109);   
rgb(490pt)=(0.176483,0.716935,0.64468);   
rgb(491pt)=(0.178806,0.717395,0.64325);   
rgb(492pt)=(0.181129,0.717854,0.641821);  
rgb(493pt)=(0.183452,0.718314,0.640391);  
rgb(494pt)=(0.185755,0.718783,0.638952);  
rgb(495pt)=(0.188027,0.719268,0.637497);  
rgb(496pt)=(0.190299,0.719753,0.636042);  
rgb(497pt)=(0.192571,0.720238,0.634587);  
rgb(498pt)=(0.194913,0.720711,0.633132);  
rgb(499pt)=(0.197338,0.72117,0.631677);   
rgb(500pt)=(0.199762,0.72163,0.630223);   
rgb(501pt)=(0.202187,0.722089,0.628768);  
rgb(502pt)=(0.204612,0.722549,0.627299);  
rgb(503pt)=(0.207037,0.723008,0.625818);  
rgb(504pt)=(0.209462,0.723468,0.624338);  
rgb(505pt)=(0.211887,0.723927,0.622857);  
rgb(506pt)=(0.214328,0.724386,0.621377);  
rgb(507pt)=(0.216778,0.724846,0.619896);  
rgb(508pt)=(0.219229,0.725305,0.618416);  
rgb(509pt)=(0.221679,0.725765,0.616935);  
rgb(510pt)=(0.224202,0.726188,0.615455);  
rgb(511pt)=(0.226754,0.726597,0.613974);  
rgb(512pt)=(0.229307,0.727005,0.612494);  
rgb(513pt)=(0.231859,0.727414,0.611014);  
rgb(514pt)=(0.234392,0.727842,0.609513);  
rgb(515pt)=(0.236919,0.728276,0.608007);  
rgb(516pt)=(0.239446,0.72871,0.606501);   
rgb(517pt)=(0.241973,0.729144,0.604995);  
rgb(518pt)=(0.244611,0.729556,0.603467);  
rgb(519pt)=(0.247266,0.729964,0.601935);  
rgb(520pt)=(0.24992,0.730372,0.600404);   
rgb(521pt)=(0.252575,0.730781,0.598872);  
rgb(522pt)=(0.25523,0.731189,0.597365);   
rgb(523pt)=(0.257884,0.731598,0.595859);  
rgb(524pt)=(0.260539,0.732006,0.594353);  
rgb(525pt)=(0.263194,0.732414,0.592846);  
rgb(526pt)=(0.265848,0.732796,0.591314);  
rgb(527pt)=(0.268503,0.733179,0.589783);  
rgb(528pt)=(0.271158,0.733562,0.588251);  
rgb(529pt)=(0.27383,0.733945,0.58672);    
rgb(530pt)=(0.276638,0.734328,0.585188);  
rgb(531pt)=(0.279446,0.734711,0.583657);  
rgb(532pt)=(0.282254,0.735094,0.582125);  
rgb(533pt)=(0.285051,0.735471,0.580594);  
rgb(534pt)=(0.287808,0.735829,0.579062);  
rgb(535pt)=(0.290565,0.736186,0.577531);  
rgb(536pt)=(0.293322,0.736544,0.575999);  
rgb(537pt)=(0.2961,0.736894,0.574468);    
rgb(538pt)=(0.298933,0.737226,0.572936);  
rgb(539pt)=(0.301767,0.737557,0.571405);  
rgb(540pt)=(0.3046,0.737889,0.569873);    
rgb(541pt)=(0.307452,0.738221,0.568351);  
rgb(542pt)=(0.310336,0.738553,0.566845);  
rgb(543pt)=(0.313221,0.738885,0.565339);  
rgb(544pt)=(0.316105,0.739217,0.563833);  
rgb(545pt)=(0.318978,0.739537,0.562315);  
rgb(546pt)=(0.321837,0.739843,0.560784);  
rgb(547pt)=(0.324696,0.74015,0.559252);   
rgb(548pt)=(0.327555,0.740456,0.557721);  
rgb(549pt)=(0.330468,0.740749,0.556216);  
rgb(550pt)=(0.333429,0.741029,0.554736);  
rgb(551pt)=(0.336389,0.74131,0.553255);   
rgb(552pt)=(0.33935,0.741591,0.551775);   
rgb(553pt)=(0.342296,0.741872,0.550279);  
rgb(554pt)=(0.345231,0.742153,0.548773);  
rgb(555pt)=(0.348167,0.742433,0.547267);  
rgb(556pt)=(0.351102,0.742714,0.545761);  
rgb(557pt)=(0.354038,0.742977,0.54429);   
rgb(558pt)=(0.356973,0.743232,0.542835);  
rgb(559pt)=(0.359908,0.743488,0.54138);   
rgb(560pt)=(0.362844,0.743743,0.539925);  
rgb(561pt)=(0.365839,0.743959,0.53847);   
rgb(562pt)=(0.368851,0.744163,0.537015);  
rgb(563pt)=(0.371863,0.744367,0.53556);   
rgb(564pt)=(0.374875,0.744571,0.534105);  
rgb(565pt)=(0.377843,0.744775,0.532672);  
rgb(566pt)=(0.380804,0.74498,0.531243);   
rgb(567pt)=(0.383765,0.745184,0.529814);  
rgb(568pt)=(0.386726,0.745388,0.528384);  
rgb(569pt)=(0.389711,0.745568,0.527003);  
rgb(570pt)=(0.392697,0.745747,0.525624);  
rgb(571pt)=(0.395684,0.745926,0.524246);  
rgb(572pt)=(0.39867,0.746104,0.522868);   
rgb(573pt)=(0.401657,0.746257,0.521489);  
rgb(574pt)=(0.404643,0.74641,0.520111);   
rgb(575pt)=(0.40763,0.746563,0.518732);   
rgb(576pt)=(0.410611,0.746716,0.517359);  
rgb(577pt)=(0.413546,0.746869,0.516032);  
rgb(578pt)=(0.416482,0.747023,0.514705);  
rgb(579pt)=(0.419417,0.747176,0.513377);  
rgb(580pt)=(0.422357,0.747319,0.512055);  
rgb(581pt)=(0.425318,0.747421,0.510753);  
rgb(582pt)=(0.428279,0.747523,0.509451);  
rgb(583pt)=(0.43124,0.747626,0.50815);    
rgb(584pt)=(0.43418,0.747735,0.506848);   
rgb(585pt)=(0.437065,0.747862,0.505546);  
rgb(586pt)=(0.439949,0.74799,0.504244);   
rgb(587pt)=(0.442834,0.748117,0.502942);  
rgb(588pt)=(0.445727,0.748227,0.501659);  
rgb(589pt)=(0.448637,0.748304,0.500408);  
rgb(590pt)=(0.451547,0.74838,0.499157);   
rgb(591pt)=(0.454457,0.748457,0.497906);  
rgb(592pt)=(0.457333,0.748522,0.496667);  
rgb(593pt)=(0.460167,0.748573,0.495441);  
rgb(594pt)=(0.463,0.748624,0.494216);     
rgb(595pt)=(0.465833,0.748675,0.492991);  
rgb(596pt)=(0.468667,0.748726,0.491779);  
rgb(597pt)=(0.4715,0.748777,0.490579);    
rgb(598pt)=(0.474333,0.748829,0.48938);   
rgb(599pt)=(0.477167,0.74888,0.48818);    
rgb(600pt)=(0.479969,0.748931,0.486995);  
rgb(601pt)=(0.482752,0.748982,0.485821);  
rgb(602pt)=(0.485534,0.749033,0.484647);  
rgb(603pt)=(0.488316,0.749084,0.483473);  
rgb(604pt)=(0.491081,0.7491,0.482316);    
rgb(605pt)=(0.493838,0.7491,0.481168);    
rgb(606pt)=(0.496595,0.7491,0.480019);    
rgb(607pt)=(0.499351,0.7491,0.47887);     
rgb(608pt)=(0.502069,0.74912,0.477722);   
rgb(609pt)=(0.504775,0.749145,0.476573);  
rgb(610pt)=(0.50748,0.749171,0.475424);   
rgb(611pt)=(0.510186,0.749196,0.474276);  
rgb(612pt)=(0.512892,0.7492,0.47317);     
rgb(613pt)=(0.515598,0.7492,0.472073);    
rgb(614pt)=(0.518303,0.7492,0.470975);    
rgb(615pt)=(0.521009,0.7492,0.469877);    
rgb(616pt)=(0.523644,0.749176,0.46878);   
rgb(617pt)=(0.526273,0.749151,0.467682);  
rgb(618pt)=(0.528902,0.749125,0.466585);  
rgb(619pt)=(0.531531,0.7491,0.465487);    
rgb(620pt)=(0.53416,0.749074,0.464415);   
rgb(621pt)=(0.536789,0.749049,0.463343);  
rgb(622pt)=(0.539418,0.749023,0.462271);  
rgb(623pt)=(0.542043,0.748998,0.461199);  
rgb(624pt)=(0.544621,0.748972,0.460127);  
rgb(625pt)=(0.547199,0.748947,0.459055);  
rgb(626pt)=(0.549777,0.748921,0.457983);  
rgb(627pt)=(0.55235,0.748891,0.45692);    
rgb(628pt)=(0.554903,0.74884,0.455899);   
rgb(629pt)=(0.557456,0.748789,0.454878);  
rgb(630pt)=(0.560008,0.748738,0.453857);  
rgb(631pt)=(0.562554,0.748687,0.452829);  
rgb(632pt)=(0.565081,0.748636,0.451783);  
rgb(633pt)=(0.567608,0.748585,0.450736);  
rgb(634pt)=(0.570135,0.748534,0.449689);  
rgb(635pt)=(0.572653,0.748474,0.44866);   
rgb(636pt)=(0.575155,0.748397,0.447665);  
rgb(637pt)=(0.577656,0.748321,0.446669);  
rgb(638pt)=(0.580158,0.748244,0.445674);  
rgb(639pt)=(0.582649,0.748168,0.444678);  
rgb(640pt)=(0.585125,0.748091,0.443683);  
rgb(641pt)=(0.587601,0.748014,0.442687);  
rgb(642pt)=(0.590077,0.747938,0.441692);  
rgb(643pt)=(0.59254,0.747861,0.440709);   
rgb(644pt)=(0.59499,0.747785,0.439739);   
rgb(645pt)=(0.597441,0.747708,0.438769);  
rgb(646pt)=(0.599891,0.747632,0.437799);  
rgb(647pt)=(0.602311,0.747555,0.436814);  
rgb(648pt)=(0.604711,0.747478,0.435819);  
rgb(649pt)=(0.60711,0.747402,0.434823);   
rgb(650pt)=(0.60951,0.747325,0.433828);   
rgb(651pt)=(0.611909,0.747232,0.432867);  
rgb(652pt)=(0.614308,0.747129,0.431922);  
rgb(653pt)=(0.616708,0.747027,0.430978);  
rgb(654pt)=(0.619107,0.746925,0.430033);  
rgb(655pt)=(0.621487,0.746823,0.429089);  
rgb(656pt)=(0.623861,0.746721,0.428144);  
rgb(657pt)=(0.626235,0.746619,0.4272);    
rgb(658pt)=(0.628609,0.746517,0.426256);  
rgb(659pt)=(0.630962,0.746393,0.425311);  
rgb(660pt)=(0.63331,0.746266,0.424367);   
rgb(661pt)=(0.635658,0.746138,0.423422);  
rgb(662pt)=(0.638007,0.746011,0.422478);  
rgb(663pt)=(0.640332,0.745906,0.421557);  
rgb(664pt)=(0.642654,0.745804,0.420638);  
rgb(665pt)=(0.644977,0.745702,0.419719);  
rgb(666pt)=(0.6473,0.7456,0.4188);        
rgb(667pt)=(0.649623,0.745472,0.417881);  
rgb(668pt)=(0.651946,0.745345,0.416962);  
rgb(669pt)=(0.654268,0.745217,0.416043);  
rgb(670pt)=(0.656587,0.745089,0.415124);  
rgb(671pt)=(0.658859,0.744962,0.414205);  
rgb(672pt)=(0.661131,0.744834,0.413286);  
rgb(673pt)=(0.663402,0.744707,0.412368);  
rgb(674pt)=(0.665674,0.744579,0.411453);  
rgb(675pt)=(0.667946,0.744451,0.410559);  
rgb(676pt)=(0.670218,0.744324,0.409666);  
rgb(677pt)=(0.672489,0.744196,0.408773);  
rgb(678pt)=(0.674755,0.744062,0.407879);  
rgb(679pt)=(0.677001,0.743909,0.406986);  
rgb(680pt)=(0.679247,0.743756,0.406092);  
rgb(681pt)=(0.681494,0.743603,0.405199);  
rgb(682pt)=(0.68374,0.743458,0.404306);   
rgb(683pt)=(0.685986,0.74333,0.403412);   
rgb(684pt)=(0.688232,0.743203,0.402519);  
rgb(685pt)=(0.690479,0.743075,0.401626);  
rgb(686pt)=(0.692704,0.742937,0.400732);  
rgb(687pt)=(0.694899,0.742784,0.399839);  
rgb(688pt)=(0.697094,0.742631,0.398945);  
rgb(689pt)=(0.699289,0.742477,0.398052);  
rgb(690pt)=(0.701497,0.742324,0.397171);  
rgb(691pt)=(0.703718,0.742171,0.396303);  
rgb(692pt)=(0.705939,0.742018,0.395435);  
rgb(693pt)=(0.708159,0.741865,0.394568);  
rgb(694pt)=(0.710351,0.741712,0.3937);    
rgb(695pt)=(0.71252,0.741559,0.392832);   
rgb(696pt)=(0.71469,0.741405,0.391964);   
rgb(697pt)=(0.71686,0.741252,0.391096);   
rgb(698pt)=(0.719029,0.741082,0.390228);  
rgb(699pt)=(0.721199,0.740904,0.38936);   
rgb(700pt)=(0.723369,0.740725,0.388492);  
rgb(701pt)=(0.725538,0.740546,0.387625);  
rgb(702pt)=(0.727708,0.740386,0.386757);  
rgb(703pt)=(0.729878,0.740233,0.385889);  
rgb(704pt)=(0.732047,0.74008,0.385021);   
rgb(705pt)=(0.734217,0.739927,0.384153);  
rgb(706pt)=(0.736366,0.739753,0.383285);  
rgb(707pt)=(0.73851,0.739574,0.382417);   
rgb(708pt)=(0.740654,0.739395,0.38155);   
rgb(709pt)=(0.742798,0.739217,0.380682);  
rgb(710pt)=(0.744919,0.739038,0.379837);  
rgb(711pt)=(0.747038,0.738859,0.378995);  
rgb(712pt)=(0.749156,0.738681,0.378152);  
rgb(713pt)=(0.751275,0.738502,0.37731);   
rgb(714pt)=(0.753394,0.738323,0.376442);  
rgb(715pt)=(0.755512,0.738145,0.375574);  
rgb(716pt)=(0.757631,0.737966,0.374707);  
rgb(717pt)=(0.75975,0.737789,0.373841);   
rgb(718pt)=(0.761868,0.737636,0.372998);  
rgb(719pt)=(0.763987,0.737483,0.372156);  
rgb(720pt)=(0.766105,0.73733,0.371314);   
rgb(721pt)=(0.76822,0.737169,0.370471);   
rgb(722pt)=(0.770313,0.736965,0.369629);  
rgb(723pt)=(0.772406,0.73676,0.368786);   
rgb(724pt)=(0.774499,0.736556,0.367944);  
rgb(725pt)=(0.776592,0.736358,0.367096);  
rgb(726pt)=(0.778686,0.736179,0.366228);  
rgb(727pt)=(0.780779,0.736001,0.36536);   
rgb(728pt)=(0.782872,0.735822,0.364492);  
rgb(729pt)=(0.784957,0.735643,0.363632);  
rgb(730pt)=(0.787024,0.735465,0.36279);   
rgb(731pt)=(0.789092,0.735286,0.361948);  
rgb(732pt)=(0.791159,0.735107,0.361105);  
rgb(733pt)=(0.793227,0.734929,0.360263);  
rgb(734pt)=(0.795295,0.73475,0.359421);   
rgb(735pt)=(0.797362,0.734571,0.358578);  
rgb(736pt)=(0.79943,0.734392,0.357736);   
rgb(737pt)=(0.801485,0.734214,0.356881);  
rgb(738pt)=(0.803527,0.734035,0.356014);  
rgb(739pt)=(0.805569,0.733856,0.355146);  
rgb(740pt)=(0.807611,0.733678,0.354278);  
rgb(741pt)=(0.809668,0.733499,0.353424);  
rgb(742pt)=(0.811735,0.73332,0.352582);   
rgb(743pt)=(0.813803,0.733142,0.35174);   
rgb(744pt)=(0.81587,0.732963,0.350897);   
rgb(745pt)=(0.817921,0.732784,0.350038);  
rgb(746pt)=(0.819963,0.732606,0.349171);  
rgb(747pt)=(0.822005,0.732427,0.348303);  
rgb(748pt)=(0.824047,0.732248,0.347435);  
rgb(749pt)=(0.826071,0.73207,0.346567);   
rgb(750pt)=(0.828087,0.731891,0.345699);  
rgb(751pt)=(0.830104,0.731712,0.344831);  
rgb(752pt)=(0.83212,0.731534,0.343963);   
rgb(753pt)=(0.834158,0.731355,0.343095);  
rgb(754pt)=(0.8362,0.731176,0.342228);    
rgb(755pt)=(0.838242,0.730998,0.34136);   
rgb(756pt)=(0.840284,0.730819,0.340492);  
rgb(757pt)=(0.842303,0.73064,0.339624);   
rgb(758pt)=(0.84432,0.730462,0.338756);   
rgb(759pt)=(0.846336,0.730283,0.337888);  
rgb(760pt)=(0.848353,0.730104,0.33702);   
rgb(761pt)=(0.850369,0.729926,0.336153);  
rgb(762pt)=(0.852386,0.729747,0.335285);  
rgb(763pt)=(0.854402,0.729568,0.334417);  
rgb(764pt)=(0.856419,0.729391,0.333546);  
rgb(765pt)=(0.858435,0.729238,0.332627);  
rgb(766pt)=(0.860452,0.729085,0.331708);  
rgb(767pt)=(0.862468,0.728932,0.330789);  
rgb(768pt)=(0.864481,0.728778,0.329874);  
rgb(769pt)=(0.866472,0.728625,0.32898);   
rgb(770pt)=(0.868463,0.728472,0.328087);  
rgb(771pt)=(0.870454,0.728319,0.327194);  
rgb(772pt)=(0.872445,0.728166,0.326295);  
rgb(773pt)=(0.874436,0.728013,0.325376);  
rgb(774pt)=(0.876427,0.727859,0.324457);  
rgb(775pt)=(0.878418,0.727706,0.323538);  
rgb(776pt)=(0.880417,0.727561,0.322619);  
rgb(777pt)=(0.882433,0.727433,0.3217);    
rgb(778pt)=(0.88445,0.727306,0.320781);   
rgb(779pt)=(0.886466,0.727178,0.319862);  
rgb(780pt)=(0.888463,0.72705,0.318933);   
rgb(781pt)=(0.890429,0.726923,0.317989);  
rgb(782pt)=(0.892394,0.726795,0.317044);  
rgb(783pt)=(0.894359,0.726668,0.3161);    
rgb(784pt)=(0.896337,0.726552,0.315132);  
rgb(785pt)=(0.898328,0.72645,0.314136);   
rgb(786pt)=(0.900319,0.726348,0.313141);  
rgb(787pt)=(0.90231,0.726246,0.312145);   
rgb(788pt)=(0.904301,0.726158,0.31115);   
rgb(789pt)=(0.906292,0.726081,0.310154);  
rgb(790pt)=(0.908283,0.726005,0.309159);  
rgb(791pt)=(0.910274,0.725928,0.308163);  
rgb(792pt)=(0.912249,0.725851,0.307151);  
rgb(793pt)=(0.914214,0.725775,0.30613);   
rgb(794pt)=(0.91618,0.725698,0.305109);   
rgb(795pt)=(0.918145,0.725622,0.304088);  
rgb(796pt)=(0.920111,0.7256,0.303031);    
rgb(797pt)=(0.922076,0.7256,0.301959);    
rgb(798pt)=(0.924041,0.7256,0.300886);    
rgb(799pt)=(0.926007,0.7256,0.299814);    
rgb(800pt)=(0.927972,0.7256,0.298722);    
rgb(801pt)=(0.929938,0.7256,0.297624);    
rgb(802pt)=(0.931903,0.7256,0.296527);    
rgb(803pt)=(0.933869,0.7256,0.295429);    
rgb(804pt)=(0.935812,0.725668,0.294264);  
rgb(805pt)=(0.937752,0.725744,0.29309);   
rgb(806pt)=(0.939692,0.725821,0.291916);  
rgb(807pt)=(0.941632,0.725897,0.290741);  
rgb(808pt)=(0.943571,0.726023,0.289518);  
rgb(809pt)=(0.945511,0.726151,0.288293);  
rgb(810pt)=(0.947451,0.726278,0.287068);  
rgb(811pt)=(0.949389,0.726411,0.285839);  
rgb(812pt)=(0.951278,0.726641,0.284537);  
rgb(813pt)=(0.953167,0.72687,0.283235);   
rgb(814pt)=(0.955056,0.7271,0.281933);    
rgb(815pt)=(0.956938,0.72734,0.280622);   
rgb(816pt)=(0.958776,0.727646,0.279243);  
rgb(817pt)=(0.960614,0.727952,0.277865);  
rgb(818pt)=(0.962451,0.728259,0.276486);  
rgb(819pt)=(0.964273,0.728597,0.275086);  
rgb(820pt)=(0.966034,0.729057,0.273606);  
rgb(821pt)=(0.967795,0.729516,0.272126);  
rgb(822pt)=(0.969557,0.729976,0.270645);  
rgb(823pt)=(0.971288,0.730473,0.269135);  
rgb(824pt)=(0.972947,0.73106,0.267552);   
rgb(825pt)=(0.974606,0.731647,0.265969);  
rgb(826pt)=(0.976265,0.732234,0.264387);  
rgb(827pt)=(0.977857,0.732879,0.262785);  
rgb(828pt)=(0.979338,0.733619,0.261151);  
rgb(829pt)=(0.980818,0.734359,0.259518);  
rgb(830pt)=(0.982299,0.735099,0.257884);  
rgb(831pt)=(0.983697,0.73591,0.256227);   
rgb(832pt)=(0.984999,0.736803,0.254542);  
rgb(833pt)=(0.986301,0.737697,0.252858);  
rgb(834pt)=(0.987603,0.73859,0.251173);   
rgb(835pt)=(0.988753,0.739566,0.249474);  
rgb(836pt)=(0.989774,0.740613,0.247764);  
rgb(837pt)=(0.990795,0.741659,0.246054);  
rgb(838pt)=(0.991816,0.742706,0.244344);  
rgb(839pt)=(0.992677,0.743816,0.242681);  
rgb(840pt)=(0.993443,0.744965,0.241048);  
rgb(841pt)=(0.994209,0.746114,0.239414);  
rgb(842pt)=(0.994975,0.747262,0.23778);   
rgb(843pt)=(0.995578,0.748465,0.236165);  
rgb(844pt)=(0.996114,0.74969,0.234557);   
rgb(845pt)=(0.99665,0.750915,0.232949);   
rgb(846pt)=(0.997186,0.752141,0.231341);  
rgb(847pt)=(0.997562,0.753386,0.229813);  
rgb(848pt)=(0.997893,0.754637,0.228307);  
rgb(849pt)=(0.998225,0.755887,0.226801);  
rgb(850pt)=(0.998557,0.757138,0.225295);  
rgb(851pt)=(0.998711,0.758433,0.223856);  
rgb(852pt)=(0.998839,0.759735,0.222426);  
rgb(853pt)=(0.998966,0.761037,0.220997);  
rgb(854pt)=(0.999094,0.762339,0.219567);  
rgb(855pt)=(0.999076,0.763641,0.218186);  
rgb(856pt)=(0.99905,0.764942,0.216808);   
rgb(857pt)=(0.999025,0.766244,0.21543);   
rgb(858pt)=(0.998995,0.767546,0.214054);  
rgb(859pt)=(0.998868,0.768848,0.212752);  
rgb(860pt)=(0.99874,0.77015,0.21145);     
rgb(861pt)=(0.998613,0.771451,0.210149);  
rgb(862pt)=(0.998473,0.772756,0.208856);  
rgb(863pt)=(0.998243,0.774083,0.207631);  
rgb(864pt)=(0.998014,0.775411,0.206405);  
rgb(865pt)=(0.997784,0.776738,0.20518);   
rgb(866pt)=(0.997539,0.77806,0.20396);    
rgb(867pt)=(0.997232,0.779362,0.20276);   
rgb(868pt)=(0.996926,0.780664,0.201561);  
rgb(869pt)=(0.99662,0.781966,0.200361);   
rgb(870pt)=(0.996299,0.783268,0.199168);  
rgb(871pt)=(0.995942,0.784569,0.197994);  
rgb(872pt)=(0.995584,0.785871,0.19682);   
rgb(873pt)=(0.995227,0.787173,0.195646);  
rgb(874pt)=(0.994842,0.788475,0.19449);   
rgb(875pt)=(0.994408,0.789777,0.193367);  
rgb(876pt)=(0.993974,0.791078,0.192244);  
rgb(877pt)=(0.99354,0.79238,0.191121);    
rgb(878pt)=(0.993083,0.793671,0.190021);  
rgb(879pt)=(0.992598,0.794947,0.188949);  
rgb(880pt)=(0.992113,0.796223,0.187877);  
rgb(881pt)=(0.991628,0.797499,0.186805);  
rgb(882pt)=(0.99113,0.798789,0.185732);   
rgb(883pt)=(0.990619,0.800091,0.18466);   
rgb(884pt)=(0.990109,0.801393,0.183588);  
rgb(885pt)=(0.989598,0.802695,0.182516);  
rgb(886pt)=(0.989072,0.803996,0.18146);   
rgb(887pt)=(0.988536,0.805298,0.180413);  
rgb(888pt)=(0.988,0.8066,0.179367);       
rgb(889pt)=(0.987464,0.807902,0.17832);   
rgb(890pt)=(0.98691,0.809186,0.177291);   
rgb(891pt)=(0.986349,0.810462,0.17627);   
rgb(892pt)=(0.985787,0.811738,0.175249);  
rgb(893pt)=(0.985226,0.813015,0.174228);  
rgb(894pt)=(0.984644,0.814311,0.173207);  
rgb(895pt)=(0.984057,0.815613,0.172186);  
rgb(896pt)=(0.98347,0.816914,0.171165);   
rgb(897pt)=(0.982883,0.818216,0.170144);  
rgb(898pt)=(0.982296,0.819518,0.169145);  
rgb(899pt)=(0.981709,0.82082,0.16815);    
rgb(900pt)=(0.981122,0.822122,0.167154);  
rgb(901pt)=(0.980535,0.823423,0.166159);  
rgb(902pt)=(0.979947,0.824725,0.165163);  
rgb(903pt)=(0.97936,0.826027,0.164168);   
rgb(904pt)=(0.978773,0.827329,0.163172);  
rgb(905pt)=(0.978186,0.828631,0.162177);  
rgb(906pt)=(0.977599,0.829932,0.161181);  
rgb(907pt)=(0.977012,0.831234,0.160186);  
rgb(908pt)=(0.976425,0.832536,0.15919);   
rgb(909pt)=(0.975838,0.833841,0.158195);  
rgb(910pt)=(0.975251,0.835168,0.157199);  
rgb(911pt)=(0.974664,0.836495,0.156204);  
rgb(912pt)=(0.974077,0.837823,0.155208);  
rgb(913pt)=(0.973489,0.83915,0.154213);   
rgb(914pt)=(0.972902,0.840477,0.153217);  
rgb(915pt)=(0.972315,0.841805,0.152222);  
rgb(916pt)=(0.971728,0.843132,0.151226);  
rgb(917pt)=(0.971155,0.844466,0.150224);  
rgb(918pt)=(0.970619,0.845819,0.149203);  
rgb(919pt)=(0.970083,0.847172,0.148182);  
rgb(920pt)=(0.969547,0.848525,0.147161);  
rgb(921pt)=(0.96902,0.849886,0.14614);    
rgb(922pt)=(0.968509,0.851265,0.145119);  
rgb(923pt)=(0.967999,0.852643,0.144098);  
rgb(924pt)=(0.967488,0.854022,0.143077);  
rgb(925pt)=(0.967,0.855411,0.142056);     
rgb(926pt)=(0.966541,0.856815,0.141035);  
rgb(927pt)=(0.966081,0.858219,0.140014);  
rgb(928pt)=(0.965622,0.859623,0.138992);  
rgb(929pt)=(0.965189,0.86104,0.137945);   
rgb(930pt)=(0.96478,0.862469,0.136873);   
rgb(931pt)=(0.964372,0.863899,0.135801);  
rgb(932pt)=(0.963963,0.865328,0.134729);  
rgb(933pt)=(0.96357,0.866773,0.133657);   
rgb(934pt)=(0.963187,0.868228,0.132585);  
rgb(935pt)=(0.962805,0.869683,0.131513);  
rgb(936pt)=(0.962422,0.871138,0.130441);  
rgb(937pt)=(0.962091,0.87261,0.129351);   
rgb(938pt)=(0.961785,0.874091,0.128253);  
rgb(939pt)=(0.961478,0.875571,0.127156);  
rgb(940pt)=(0.961172,0.877052,0.126058);  
rgb(941pt)=(0.960885,0.878571,0.124961);  
rgb(942pt)=(0.960605,0.880103,0.123863);  
rgb(943pt)=(0.960324,0.881634,0.122765);  
rgb(944pt)=(0.960043,0.883166,0.121668);  
rgb(945pt)=(0.959849,0.884719,0.120549);  
rgb(946pt)=(0.95967,0.886276,0.119426);   
rgb(947pt)=(0.959491,0.887833,0.118302);  
rgb(948pt)=(0.959313,0.88939,0.117179);   
rgb(949pt)=(0.959181,0.890995,0.116032);  
rgb(950pt)=(0.959054,0.892603,0.114884);  
rgb(951pt)=(0.958926,0.894211,0.113735);  
rgb(952pt)=(0.958799,0.895819,0.112587);  
rgb(953pt)=(0.958748,0.897453,0.111464);  
rgb(954pt)=(0.958697,0.899086,0.110341);  
rgb(955pt)=(0.958646,0.90072,0.109217);   
rgb(956pt)=(0.958602,0.902359,0.108089);  
rgb(957pt)=(0.958628,0.904043,0.106915);  
rgb(958pt)=(0.958653,0.905728,0.105741);  
rgb(959pt)=(0.958679,0.907413,0.104567);  
rgb(960pt)=(0.958718,0.909102,0.103393);  
rgb(961pt)=(0.95882,0.910812,0.102219);   
rgb(962pt)=(0.958922,0.912522,0.101044);  
rgb(963pt)=(0.959024,0.914232,0.0998703); 
rgb(964pt)=(0.959153,0.915962,0.0986961); 
rgb(965pt)=(0.959357,0.917749,0.0975219); 
rgb(966pt)=(0.959561,0.919536,0.0963477); 
rgb(967pt)=(0.959765,0.921323,0.0951736); 
rgb(968pt)=(0.959996,0.923118,0.0939907); 
rgb(969pt)=(0.960277,0.924931,0.092791);  
rgb(970pt)=(0.960557,0.926743,0.0915913); 
rgb(971pt)=(0.960838,0.928555,0.0903916); 
rgb(972pt)=(0.961151,0.930378,0.0891919); 
rgb(973pt)=(0.961509,0.932216,0.0879922); 
rgb(974pt)=(0.961866,0.934054,0.0867925); 
rgb(975pt)=(0.962223,0.935892,0.0855928); 
rgb(976pt)=(0.96262,0.937768,0.0843802);  
rgb(977pt)=(0.963053,0.939683,0.083155);  
rgb(978pt)=(0.963487,0.941597,0.0819297); 
rgb(979pt)=(0.963921,0.943512,0.0807045); 
rgb(980pt)=(0.964415,0.945426,0.0794643); 
rgb(981pt)=(0.964951,0.947341,0.0782135); 
rgb(982pt)=(0.965487,0.949255,0.0769628); 
rgb(983pt)=(0.966023,0.951169,0.075712);  
rgb(984pt)=(0.966594,0.953118,0.0744441); 
rgb(985pt)=(0.967181,0.955083,0.0731679); 
rgb(986pt)=(0.967768,0.957049,0.0718916); 
rgb(987pt)=(0.968355,0.959014,0.0706153); 
rgb(988pt)=(0.96898,0.960999,0.0693006);  
rgb(989pt)=(0.969619,0.96299,0.0679733);  
rgb(990pt)=(0.970257,0.964981,0.0666459); 
rgb(991pt)=(0.970895,0.966972,0.0653186); 
rgb(992pt)=(0.971554,0.968984,0.0639486); 
rgb(993pt)=(0.972218,0.971001,0.0625703); 
rgb(994pt)=(0.972882,0.973017,0.0611919); 
rgb(995pt)=(0.973545,0.975034,0.0598135); 
rgb(996pt)=(0.974232,0.97705,0.058318);   
rgb(997pt)=(0.974922,0.979067,0.056812);  
rgb(998pt)=(0.975611,0.981083,0.055306);  
rgb(999pt)=(0.9763,0.9831,0.0538)
}         

%\def\mycolormap{rgb(0pt)=(0.230033,0.298999,0.754002);
%                rgb(1pt)=(0.231184,0.300753,0.755593);
%                rgb(2pt)=(0.232336,0.302506,0.75718);
%                rgb(3pt)=(0.233489,0.304258,0.758764);
%                rgb(4pt)=(0.234643,0.306009,0.760343);
%                rgb(5pt)=(0.235798,0.30776,0.761918);
%                rgb(6pt)=(0.236954,0.309509,0.76349);
%                rgb(7pt)=(0.238111,0.311258,0.765057);
%                rgb(8pt)=(0.239269,0.313006,0.76662);
%                rgb(9pt)=(0.240429,0.314754,0.76818);
%                rgb(10pt)=(0.241589,0.3165,0.769735);
%                rgb(11pt)=(0.242751,0.318246,0.771286);
%                rgb(12pt)=(0.243913,0.319991,0.772833);
%                rgb(13pt)=(0.245077,0.321735,0.774376);
%                rgb(14pt)=(0.246242,0.323478,0.775915);
%                rgb(15pt)=(0.247408,0.32522,0.77745);
%                rgb(16pt)=(0.248575,0.326962,0.77898);
%                rgb(17pt)=(0.249743,0.328702,0.780507);
%                rgb(18pt)=(0.250913,0.330442,0.782029);
%                rgb(19pt)=(0.252084,0.332181,0.783547);
%                rgb(20pt)=(0.253255,0.333919,0.785061);
%                rgb(21pt)=(0.254428,0.335656,0.78657);
%                rgb(22pt)=(0.255602,0.337393,0.788075);
%                rgb(23pt)=(0.256778,0.339128,0.789576);
%                rgb(24pt)=(0.257954,0.340863,0.791073);
%                rgb(25pt)=(0.259132,0.342597,0.792565);
%                rgb(26pt)=(0.260311,0.344329,0.794054);
%                rgb(27pt)=(0.261491,0.346061,0.795537);
%                rgb(28pt)=(0.262672,0.347792,0.797017);
%                rgb(29pt)=(0.263855,0.349522,0.798492);
%                rgb(30pt)=(0.265039,0.351251,0.799963);
%                rgb(31pt)=(0.266224,0.352979,0.801429);
%                rgb(32pt)=(0.26741,0.354706,0.802891);
%                rgb(33pt)=(0.268598,0.356433,0.804348);
%                rgb(34pt)=(0.269787,0.358158,0.805801);
%                rgb(35pt)=(0.270977,0.359882,0.80725);
%                rgb(36pt)=(0.272168,0.361605,0.808694);
%                rgb(37pt)=(0.273361,0.363328,0.810134);
%                rgb(38pt)=(0.274555,0.365049,0.811569);
%                rgb(39pt)=(0.27575,0.366769,0.813);
%                rgb(40pt)=(0.276946,0.368489,0.814426);
%                rgb(41pt)=(0.278144,0.370207,0.815847);
%                rgb(42pt)=(0.279343,0.371924,0.817264);
%                rgb(43pt)=(0.280543,0.37364,0.818677);
%                rgb(44pt)=(0.281745,0.375355,0.820085);
%                rgb(45pt)=(0.282948,0.37707,0.821488);
%                rgb(46pt)=(0.284152,0.378783,0.822887);
%                rgb(47pt)=(0.285357,0.380494,0.824281);
%                rgb(48pt)=(0.286564,0.382205,0.82567);
%                rgb(49pt)=(0.287772,0.383915,0.827055);
%                rgb(50pt)=(0.288982,0.385624,0.828435);
%                rgb(51pt)=(0.290193,0.387331,0.829811);
%                rgb(52pt)=(0.291405,0.389037,0.831181);
%                rgb(53pt)=(0.292618,0.390743,0.832547);
%                rgb(54pt)=(0.293833,0.392447,0.833909);
%                rgb(55pt)=(0.295049,0.39415,0.835265);
%                rgb(56pt)=(0.296266,0.395851,0.836617);
%                rgb(57pt)=(0.297485,0.397552,0.837964);
%                rgb(58pt)=(0.298705,0.399251,0.839306);
%                rgb(59pt)=(0.299926,0.40095,0.840644);
%                rgb(60pt)=(0.301149,0.402647,0.841977);
%                rgb(61pt)=(0.302373,0.404342,0.843304);
%                rgb(62pt)=(0.303599,0.406037,0.844627);
%                rgb(63pt)=(0.304825,0.40773,0.845946);
%                rgb(64pt)=(0.306053,0.409422,0.847259);
%                rgb(65pt)=(0.307283,0.411113,0.848567);
%                rgb(66pt)=(0.308514,0.412802,0.849871);
%                rgb(67pt)=(0.309746,0.414491,0.85117);
%                rgb(68pt)=(0.310979,0.416177,0.852463);
%                rgb(69pt)=(0.312214,0.417863,0.853752);
%                rgb(70pt)=(0.31345,0.419547,0.855036);
%                rgb(71pt)=(0.314687,0.42123,0.856315);
%                rgb(72pt)=(0.315926,0.422912,0.857589);
%                rgb(73pt)=(0.317166,0.424592,0.858858);
%                rgb(74pt)=(0.318408,0.426271,0.860122);
%                rgb(75pt)=(0.319651,0.427949,0.861381);
%                rgb(76pt)=(0.320895,0.429625,0.862635);
%                rgb(77pt)=(0.32214,0.4313,0.863884);
%                rgb(78pt)=(0.323387,0.432973,0.865128);
%                rgb(79pt)=(0.324635,0.434645,0.866367);
%                rgb(80pt)=(0.325885,0.436315,0.867601);
%                rgb(81pt)=(0.327136,0.437984,0.868829);
%                rgb(82pt)=(0.328388,0.439652,0.870053);
%                rgb(83pt)=(0.329641,0.441318,0.871272);
%                rgb(84pt)=(0.330896,0.442983,0.872485);
%                rgb(85pt)=(0.332152,0.444646,0.873693);
%                rgb(86pt)=(0.33341,0.446308,0.874897);
%                rgb(87pt)=(0.334669,0.447968,0.876095);
%                rgb(88pt)=(0.335929,0.449627,0.877287);
%                rgb(89pt)=(0.33719,0.451284,0.878475);
%                rgb(90pt)=(0.338453,0.452939,0.879658);
%                rgb(91pt)=(0.339717,0.454593,0.880835);
%                rgb(92pt)=(0.340983,0.456246,0.882007);
%                rgb(93pt)=(0.342249,0.457897,0.883174);
%                rgb(94pt)=(0.343517,0.459546,0.884335);
%                rgb(95pt)=(0.344787,0.461194,0.885492);
%                rgb(96pt)=(0.346058,0.46284,0.886643);
%                rgb(97pt)=(0.347329,0.464484,0.887789);
%                rgb(98pt)=(0.348603,0.466127,0.888929);
%                rgb(99pt)=(0.349877,0.467768,0.890065);
%                rgb(100pt)=(0.351153,0.469407,0.891195);
%                rgb(101pt)=(0.35243,0.471045,0.892319);
%                rgb(102pt)=(0.353709,0.472681,0.893439);
%                rgb(103pt)=(0.354988,0.474316,0.894553);
%                rgb(104pt)=(0.356269,0.475948,0.895661);
%                rgb(105pt)=(0.357552,0.477579,0.896764);
%                rgb(106pt)=(0.358835,0.479209,0.897862);
%                rgb(107pt)=(0.36012,0.480836,0.898955);
%                rgb(108pt)=(0.361406,0.482462,0.900042);
%                rgb(109pt)=(0.362693,0.484086,0.901124);
%                rgb(110pt)=(0.363982,0.485708,0.9022);
%                rgb(111pt)=(0.365272,0.487328,0.903271);
%                rgb(112pt)=(0.366563,0.488947,0.904337);
%                rgb(113pt)=(0.367855,0.490564,0.905397);
%                rgb(114pt)=(0.369149,0.492179,0.906451);
%                rgb(115pt)=(0.370444,0.493792,0.9075);
%                rgb(116pt)=(0.37174,0.495403,0.908544);
%                rgb(117pt)=(0.373037,0.497012,0.909582);
%                rgb(118pt)=(0.374335,0.49862,0.910615);
%                rgb(119pt)=(0.375635,0.500225,0.911642);
%                rgb(120pt)=(0.376936,0.501829,0.912664);
%                rgb(121pt)=(0.378238,0.50343,0.91368);
%                rgb(122pt)=(0.379541,0.50503,0.91469);
%                rgb(123pt)=(0.380846,0.506628,0.915695);
%                rgb(124pt)=(0.382151,0.508224,0.916695);
%                rgb(125pt)=(0.383458,0.509818,0.917689);
%                rgb(126pt)=(0.384766,0.51141,0.918677);
%                rgb(127pt)=(0.386075,0.513,0.91966);
%                rgb(128pt)=(0.387386,0.514587,0.920637);
%                rgb(129pt)=(0.388697,0.516173,0.921608);
%                rgb(130pt)=(0.39001,0.517757,0.922574);
%                rgb(131pt)=(0.391324,0.519339,0.923535);
%                rgb(132pt)=(0.392639,0.520919,0.924489);
%                rgb(133pt)=(0.393955,0.522496,0.925438);
%                rgb(134pt)=(0.395272,0.524072,0.926382);
%                rgb(135pt)=(0.396591,0.525646,0.927319);
%                rgb(136pt)=(0.39791,0.527217,0.928251);
%                rgb(137pt)=(0.399231,0.528786,0.929178);
%                rgb(138pt)=(0.400552,0.530353,0.930099);
%                rgb(139pt)=(0.401875,0.531918,0.931014);
%                rgb(140pt)=(0.403199,0.533481,0.931923);
%                rgb(141pt)=(0.404524,0.535042,0.932826);
%                rgb(142pt)=(0.40585,0.5366,0.933724);
%                rgb(143pt)=(0.407177,0.538156,0.934616);
%                rgb(144pt)=(0.408505,0.53971,0.935503);
%                rgb(145pt)=(0.409834,0.541262,0.936383);
%                rgb(146pt)=(0.411165,0.542812,0.937258);
%                rgb(147pt)=(0.412496,0.544359,0.938127);
%                rgb(148pt)=(0.413828,0.545904,0.938991);
%                rgb(149pt)=(0.415162,0.547447,0.939848);
%                rgb(150pt)=(0.416496,0.548988,0.9407);
%                rgb(151pt)=(0.417832,0.550526,0.941546);
%                rgb(152pt)=(0.419168,0.552062,0.942386);
%                rgb(153pt)=(0.420505,0.553595,0.94322);
%                rgb(154pt)=(0.421844,0.555127,0.944049);
%                rgb(155pt)=(0.423183,0.556656,0.944871);
%                rgb(156pt)=(0.424524,0.558182,0.945688);
%                rgb(157pt)=(0.425865,0.559706,0.946499);
%                rgb(158pt)=(0.427207,0.561228,0.947304);
%                rgb(159pt)=(0.42855,0.562747,0.948104);
%                rgb(160pt)=(0.429895,0.564264,0.948897);
%                rgb(161pt)=(0.43124,0.565779,0.949685);
%                rgb(162pt)=(0.432586,0.567291,0.950466);
%                rgb(163pt)=(0.433933,0.568801,0.951242);
%                rgb(164pt)=(0.435281,0.570308,0.952012);
%                rgb(165pt)=(0.436629,0.571813,0.952776);
%                rgb(166pt)=(0.437979,0.573315,0.953534);
%                rgb(167pt)=(0.43933,0.574815,0.954286);
%                rgb(168pt)=(0.440681,0.576312,0.955032);
%                rgb(169pt)=(0.442033,0.577807,0.955772);
%                rgb(170pt)=(0.443386,0.579299,0.956507);
%                rgb(171pt)=(0.44474,0.580788,0.957235);
%                rgb(172pt)=(0.446095,0.582275,0.957957);
%                rgb(173pt)=(0.447451,0.58376,0.958674);
%                rgb(174pt)=(0.448807,0.585242,0.959384);
%                rgb(175pt)=(0.450165,0.586721,0.960089);
%                rgb(176pt)=(0.451523,0.588198,0.960787);
%                rgb(177pt)=(0.452882,0.589672,0.96148);
%                rgb(178pt)=(0.454241,0.591143,0.962166);
%                rgb(179pt)=(0.455602,0.592612,0.962847);
%                rgb(180pt)=(0.456963,0.594078,0.963521);
%                rgb(181pt)=(0.458325,0.595541,0.96419);
%                rgb(182pt)=(0.459688,0.597002,0.964852);
%                rgb(183pt)=(0.461051,0.59846,0.965509);
%                rgb(184pt)=(0.462415,0.599916,0.966159);
%                rgb(185pt)=(0.46378,0.601368,0.966804);
%                rgb(186pt)=(0.465146,0.602818,0.967442);
%                rgb(187pt)=(0.466512,0.604265,0.968074);
%                rgb(188pt)=(0.467879,0.605709,0.9687);
%                rgb(189pt)=(0.469247,0.607151,0.969321);
%                rgb(190pt)=(0.470615,0.60859,0.969935);
%                rgb(191pt)=(0.471984,0.610026,0.970543);
%                rgb(192pt)=(0.473354,0.611459,0.971145);
%                rgb(193pt)=(0.474724,0.612889,0.97174);
%                rgb(194pt)=(0.476095,0.614317,0.97233);
%                rgb(195pt)=(0.477467,0.615741,0.972914);
%                rgb(196pt)=(0.478839,0.617163,0.973492);
%                rgb(197pt)=(0.480212,0.618582,0.974063);
%                rgb(198pt)=(0.481585,0.619998,0.974628);
%                rgb(199pt)=(0.482959,0.621411,0.975188);
%                rgb(200pt)=(0.484333,0.622821,0.975741);
%                rgb(201pt)=(0.485708,0.624229,0.976288);
%                rgb(202pt)=(0.487084,0.625633,0.976829);
%                rgb(203pt)=(0.48846,0.627034,0.977363);
%                rgb(204pt)=(0.489837,0.628433,0.977892);
%                rgb(205pt)=(0.491214,0.629828,0.978414);
%                rgb(206pt)=(0.492592,0.631221,0.978931);
%                rgb(207pt)=(0.49397,0.63261,0.979441);
%                rgb(208pt)=(0.495349,0.633997,0.979945);
%                rgb(209pt)=(0.496728,0.63538,0.980443);
%                rgb(210pt)=(0.498107,0.636761,0.980934);
%                rgb(211pt)=(0.499487,0.638138,0.98142);
%                rgb(212pt)=(0.500868,0.639512,0.981899);
%                rgb(213pt)=(0.502249,0.640883,0.982372);
%                rgb(214pt)=(0.50363,0.642252,0.982839);
%                rgb(215pt)=(0.505012,0.643617,0.9833);
%                rgb(216pt)=(0.506394,0.644979,0.983754);
%                rgb(217pt)=(0.507777,0.646337,0.984203);
%                rgb(218pt)=(0.509159,0.647693,0.984645);
%                rgb(219pt)=(0.510543,0.649046,0.985081);
%                rgb(220pt)=(0.511926,0.650395,0.985511);
%                rgb(221pt)=(0.51331,0.651741,0.985934);
%                rgb(222pt)=(0.514695,0.653084,0.986352);
%                rgb(223pt)=(0.516079,0.654424,0.986763);
%                rgb(224pt)=(0.517464,0.655761,0.987168);
%                rgb(225pt)=(0.518849,0.657094,0.987566);
%                rgb(226pt)=(0.520235,0.658425,0.987959);
%                rgb(227pt)=(0.521621,0.659752,0.988345);
%                rgb(228pt)=(0.523007,0.661075,0.988725);
%                rgb(229pt)=(0.524393,0.662396,0.989099);
%                rgb(230pt)=(0.525779,0.663713,0.989466);
%                rgb(231pt)=(0.527166,0.665027,0.989827);
%                rgb(232pt)=(0.528553,0.666337,0.990182);
%                rgb(233pt)=(0.52994,0.667645,0.990531);
%                rgb(234pt)=(0.531328,0.668949,0.990874);
%                rgb(235pt)=(0.532715,0.670249,0.99121);
%                rgb(236pt)=(0.534103,0.671547,0.99154);
%                rgb(237pt)=(0.535491,0.67284,0.991864);
%                rgb(238pt)=(0.536879,0.674131,0.992181);
%                rgb(239pt)=(0.538267,0.675418,0.992492);
%                rgb(240pt)=(0.539655,0.676702,0.992797);
%                rgb(241pt)=(0.541043,0.677982,0.993096);
%                rgb(242pt)=(0.542432,0.679259,0.993388);
%                rgb(243pt)=(0.54382,0.680533,0.993674);
%                rgb(244pt)=(0.545209,0.681803,0.993954);
%                rgb(245pt)=(0.546598,0.683069,0.994228);
%                rgb(246pt)=(0.547986,0.684333,0.994495);
%                rgb(247pt)=(0.549375,0.685592,0.994756);
%                rgb(248pt)=(0.550764,0.686848,0.995011);
%                rgb(249pt)=(0.552153,0.688101,0.995259);
%                rgb(250pt)=(0.553541,0.68935,0.995501);
%                rgb(251pt)=(0.55493,0.690596,0.995737);
%                rgb(252pt)=(0.556319,0.691838,0.995967);
%                rgb(253pt)=(0.557708,0.693077,0.99619);
%                rgb(254pt)=(0.559096,0.694312,0.996407);
%                rgb(255pt)=(0.560485,0.695543,0.996618);
%                rgb(256pt)=(0.561873,0.696771,0.996822);
%                rgb(257pt)=(0.563262,0.697995,0.99702);
%                rgb(258pt)=(0.56465,0.699216,0.997212);
%                rgb(259pt)=(0.566038,0.700433,0.997398);
%                rgb(260pt)=(0.567426,0.701647,0.997577);
%                rgb(261pt)=(0.568814,0.702857,0.99775);
%                rgb(262pt)=(0.570202,0.704063,0.997916);
%                rgb(263pt)=(0.57159,0.705265,0.998077);
%                rgb(264pt)=(0.572977,0.706464,0.998231);
%                rgb(265pt)=(0.574364,0.707659,0.998378);
%                rgb(266pt)=(0.575751,0.708851,0.99852);
%                rgb(267pt)=(0.577138,0.710038,0.998655);
%                rgb(268pt)=(0.578525,0.711223,0.998784);
%                rgb(269pt)=(0.579911,0.712403,0.998906);
%                rgb(270pt)=(0.581297,0.713579,0.999022);
%                rgb(271pt)=(0.582683,0.714752,0.999132);
%                rgb(272pt)=(0.584069,0.715921,0.999236);
%                rgb(273pt)=(0.585454,0.717087,0.999333);
%                rgb(274pt)=(0.586839,0.718248,0.999424);
%                rgb(275pt)=(0.588224,0.719406,0.999508);
%                rgb(276pt)=(0.589608,0.72056,0.999587);
%                rgb(277pt)=(0.590992,0.72171,0.999659);
%                rgb(278pt)=(0.592376,0.722856,0.999724);
%                rgb(279pt)=(0.593759,0.723999,0.999784);
%                rgb(280pt)=(0.595142,0.725137,0.999837);
%                rgb(281pt)=(0.596524,0.726272,0.999884);
%                rgb(282pt)=(0.597906,0.727403,0.999924);
%                rgb(283pt)=(0.599288,0.72853,0.999958);
%                rgb(284pt)=(0.600669,0.729653,0.999986);
%                rgb(285pt)=(0.602045,0.730767,1);
%                rgb(286pt)=(0.603416,0.731871,1);
%                rgb(287pt)=(0.604791,0.732976,1);
%                rgb(288pt)=(0.606169,0.734082,1);
%                rgb(289pt)=(0.60755,0.735188,1);
%                rgb(290pt)=(0.608934,0.736295,1);
%                rgb(291pt)=(0.610322,0.737402,1);
%                rgb(292pt)=(0.611702,0.738497,0.999982);
%                rgb(293pt)=(0.613078,0.739585,0.999953);
%                rgb(294pt)=(0.614455,0.740668,0.999918);
%                rgb(295pt)=(0.61583,0.741748,0.999876);
%                rgb(296pt)=(0.617205,0.742824,0.999829);
%                rgb(297pt)=(0.618579,0.743895,0.999774);
%                rgb(298pt)=(0.619953,0.744963,0.999714);
%                rgb(299pt)=(0.621326,0.746026,0.999647);
%                rgb(300pt)=(0.622699,0.747086,0.999575);
%                rgb(301pt)=(0.624071,0.748141,0.999495);
%                rgb(302pt)=(0.625442,0.749192,0.99941);
%                rgb(303pt)=(0.626812,0.75024,0.999318);
%                rgb(304pt)=(0.628182,0.751283,0.99922);
%                rgb(305pt)=(0.629551,0.752322,0.999116);
%                rgb(306pt)=(0.630919,0.753357,0.999005);
%                rgb(307pt)=(0.632287,0.754387,0.998888);
%                rgb(308pt)=(0.633654,0.755414,0.998765);
%                rgb(309pt)=(0.63502,0.756436,0.998635);
%                rgb(310pt)=(0.636385,0.757455,0.9985);
%                rgb(311pt)=(0.63775,0.758469,0.998358);
%                rgb(312pt)=(0.639113,0.759479,0.998209);
%                rgb(313pt)=(0.640476,0.760485,0.998055);
%                rgb(314pt)=(0.641838,0.761486,0.997894);
%                rgb(315pt)=(0.6432,0.762484,0.997727);
%                rgb(316pt)=(0.64456,0.763477,0.997554);
%                rgb(317pt)=(0.645919,0.764466,0.997374);
%                rgb(318pt)=(0.647278,0.76545,0.997188);
%                rgb(319pt)=(0.648636,0.766431,0.996996);
%                rgb(320pt)=(0.649992,0.767407,0.996798);
%                rgb(321pt)=(0.651348,0.768379,0.996593);
%                rgb(322pt)=(0.652703,0.769346,0.996383);
%                rgb(323pt)=(0.654057,0.77031,0.996166);
%                rgb(324pt)=(0.65541,0.771269,0.995942);
%                rgb(325pt)=(0.656762,0.772223,0.995713);
%                rgb(326pt)=(0.658113,0.773174,0.995477);
%                rgb(327pt)=(0.659463,0.77412,0.995235);
%                rgb(328pt)=(0.660812,0.775061,0.994987);
%                rgb(329pt)=(0.66216,0.775999,0.994733);
%                rgb(330pt)=(0.663507,0.776932,0.994472);
%                rgb(331pt)=(0.664853,0.77786,0.994205);
%                rgb(332pt)=(0.666197,0.778785,0.993932);
%                rgb(333pt)=(0.667541,0.779704,0.993653);
%                rgb(334pt)=(0.668884,0.78062,0.993367);
%                rgb(335pt)=(0.670225,0.781531,0.993076);
%                rgb(336pt)=(0.671565,0.782438,0.992778);
%                rgb(337pt)=(0.672904,0.78334,0.992474);
%                rgb(338pt)=(0.674242,0.784238,0.992164);
%                rgb(339pt)=(0.675579,0.785131,0.991848);
%                rgb(340pt)=(0.676915,0.78602,0.991525);
%                rgb(341pt)=(0.678249,0.786905,0.991196);
%                rgb(342pt)=(0.679583,0.787784,0.990861);
%                rgb(343pt)=(0.680915,0.78866,0.99052);
%                rgb(344pt)=(0.682245,0.789531,0.990173);
%                rgb(345pt)=(0.683575,0.790398,0.98982);
%                rgb(346pt)=(0.684903,0.79126,0.98946);
%                rgb(347pt)=(0.68623,0.792117,0.989095);
%                rgb(348pt)=(0.687556,0.79297,0.988723);
%                rgb(349pt)=(0.68888,0.793819,0.988345);
%                rgb(350pt)=(0.690203,0.794662,0.987961);
%                rgb(351pt)=(0.691525,0.795502,0.987571);
%                rgb(352pt)=(0.692845,0.796337,0.987174);
%                rgb(353pt)=(0.694164,0.797167,0.986772);
%                rgb(354pt)=(0.695481,0.797992,0.986363);
%                rgb(355pt)=(0.696798,0.798814,0.985949);
%                rgb(356pt)=(0.698112,0.79963,0.985528);
%                rgb(357pt)=(0.699426,0.800442,0.985101);
%                rgb(358pt)=(0.700738,0.801249,0.984668);
%                rgb(359pt)=(0.702048,0.802052,0.984229);
%                rgb(360pt)=(0.703357,0.80285,0.983784);
%                rgb(361pt)=(0.704665,0.803643,0.983333);
%                rgb(362pt)=(0.705971,0.804432,0.982875);
%                rgb(363pt)=(0.707275,0.805216,0.982412);
%                rgb(364pt)=(0.708578,0.805995,0.981943);
%                rgb(365pt)=(0.70988,0.80677,0.981467);
%                rgb(366pt)=(0.71118,0.80754,0.980986);
%                rgb(367pt)=(0.712478,0.808306,0.980498);
%                rgb(368pt)=(0.713775,0.809066,0.980004);
%                rgb(369pt)=(0.715071,0.809822,0.979505);
%                rgb(370pt)=(0.716364,0.810574,0.978999);
%                rgb(371pt)=(0.717656,0.81132,0.978487);
%                rgb(372pt)=(0.718947,0.812062,0.97797);
%                rgb(373pt)=(0.720235,0.812799,0.977446);
%                rgb(374pt)=(0.721523,0.813531,0.976916);
%                rgb(375pt)=(0.722808,0.814259,0.976381);
%                rgb(376pt)=(0.724092,0.814982,0.975839);
%                rgb(377pt)=(0.725374,0.8157,0.975291);
%                rgb(378pt)=(0.726654,0.816413,0.974737);
%                rgb(379pt)=(0.727933,0.817121,0.974178);
%                rgb(380pt)=(0.72921,0.817825,0.973612);
%                rgb(381pt)=(0.730485,0.818524,0.973041);
%                rgb(382pt)=(0.731759,0.819218,0.972463);
%                rgb(383pt)=(0.73303,0.819907,0.97188);
%                rgb(384pt)=(0.7343,0.820592,0.97129);
%                rgb(385pt)=(0.735568,0.821271,0.970695);
%                rgb(386pt)=(0.736834,0.821946,0.970094);
%                rgb(387pt)=(0.738099,0.822616,0.969486);
%                rgb(388pt)=(0.739361,0.823281,0.968873);
%                rgb(389pt)=(0.740622,0.823941,0.968254);
%                rgb(390pt)=(0.741881,0.824597,0.967629);
%                rgb(391pt)=(0.743138,0.825247,0.966998);
%                rgb(392pt)=(0.744393,0.825893,0.966362);
%                rgb(393pt)=(0.745646,0.826533,0.965719);
%                rgb(394pt)=(0.746897,0.827169,0.965071);
%                rgb(395pt)=(0.748147,0.8278,0.964416);
%                rgb(396pt)=(0.749394,0.828426,0.963756);
%                rgb(397pt)=(0.750639,0.829047,0.96309);
%                rgb(398pt)=(0.751883,0.829663,0.962418);
%                rgb(399pt)=(0.753124,0.830274,0.961741);
%                rgb(400pt)=(0.754363,0.83088,0.961057);
%                rgb(401pt)=(0.755601,0.831481,0.960368);
%                rgb(402pt)=(0.756836,0.832078,0.959673);
%                rgb(403pt)=(0.758069,0.832669,0.958972);
%                rgb(404pt)=(0.759301,0.833255,0.958265);
%                rgb(405pt)=(0.76053,0.833837,0.957552);
%                rgb(406pt)=(0.761757,0.834413,0.956834);
%                rgb(407pt)=(0.762982,0.834985,0.95611);
%                rgb(408pt)=(0.764205,0.835551,0.95538);
%                rgb(409pt)=(0.765425,0.836112,0.954644);
%                rgb(410pt)=(0.766644,0.836669,0.953903);
%                rgb(411pt)=(0.76786,0.83722,0.953156);
%                rgb(412pt)=(0.769075,0.837766,0.952403);
%                rgb(413pt)=(0.770287,0.838308,0.951644);
%                rgb(414pt)=(0.771497,0.838844,0.95088);
%                rgb(415pt)=(0.772704,0.839375,0.95011);
%                rgb(416pt)=(0.77391,0.839901,0.949334);
%                rgb(417pt)=(0.775113,0.840422,0.948553);
%                rgb(418pt)=(0.776314,0.840938,0.947766);
%                rgb(419pt)=(0.777512,0.841449,0.946973);
%                rgb(420pt)=(0.778709,0.841955,0.946175);
%                rgb(421pt)=(0.779903,0.842456,0.945371);
%                rgb(422pt)=(0.781095,0.842951,0.944561);
%                rgb(423pt)=(0.782284,0.843442,0.943746);
%                rgb(424pt)=(0.783471,0.843927,0.942925);
%                rgb(425pt)=(0.784656,0.844408,0.942098);
%                rgb(426pt)=(0.785839,0.844883,0.941266);
%                rgb(427pt)=(0.787019,0.845353,0.940428);
%                rgb(428pt)=(0.788196,0.845818,0.939584);
%                rgb(429pt)=(0.789372,0.846278,0.938735);
%                rgb(430pt)=(0.790545,0.846732,0.937881);
%                rgb(431pt)=(0.791715,0.847182,0.937021);
%                rgb(432pt)=(0.792883,0.847626,0.936155);
%                rgb(433pt)=(0.794049,0.848065,0.935283);
%                rgb(434pt)=(0.795212,0.848499,0.934407);
%                rgb(435pt)=(0.796373,0.848928,0.933524);
%                rgb(436pt)=(0.797531,0.849352,0.932636);
%                rgb(437pt)=(0.798686,0.84977,0.931743);
%                rgb(438pt)=(0.799839,0.850184,0.930844);
%                rgb(439pt)=(0.80099,0.850592,0.92994);
%                rgb(440pt)=(0.802138,0.850995,0.92903);
%                rgb(441pt)=(0.803284,0.851392,0.928114);
%                rgb(442pt)=(0.804427,0.851785,0.927194);
%                rgb(443pt)=(0.805567,0.852172,0.926267);
%                rgb(444pt)=(0.806705,0.852554,0.925336);
%                rgb(445pt)=(0.80784,0.852931,0.924398);
%                rgb(446pt)=(0.808972,0.853302,0.923456);
%                rgb(447pt)=(0.810102,0.853669,0.922508);
%                rgb(448pt)=(0.81123,0.85403,0.921554);
%                rgb(449pt)=(0.812354,0.854386,0.920595);
%                rgb(450pt)=(0.813476,0.854736,0.919631);
%                rgb(451pt)=(0.814595,0.855081,0.918662);
%                rgb(452pt)=(0.815712,0.855421,0.917687);
%                rgb(453pt)=(0.816826,0.855756,0.916706);
%                rgb(454pt)=(0.817937,0.856086,0.915721);
%                rgb(455pt)=(0.819045,0.85641,0.91473);
%                rgb(456pt)=(0.820151,0.856729,0.913733);
%                rgb(457pt)=(0.821254,0.857042,0.912732);
%                rgb(458pt)=(0.822354,0.85735,0.911725);
%                rgb(459pt)=(0.823451,0.857653,0.910713);
%                rgb(460pt)=(0.824546,0.857951,0.909695);
%                rgb(461pt)=(0.825637,0.858243,0.908672);
%                rgb(462pt)=(0.826726,0.85853,0.907644);
%                rgb(463pt)=(0.827812,0.858812,0.906611);
%                rgb(464pt)=(0.828895,0.859089,0.905572);
%                rgb(465pt)=(0.829975,0.85936,0.904528);
%                rgb(466pt)=(0.831053,0.859625,0.903479);
%                rgb(467pt)=(0.832127,0.859886,0.902425);
%                rgb(468pt)=(0.833199,0.860141,0.901366);
%                rgb(469pt)=(0.834268,0.86039,0.900301);
%                rgb(470pt)=(0.835333,0.860635,0.899231);
%                rgb(471pt)=(0.836396,0.860873,0.898156);
%                rgb(472pt)=(0.837456,0.861107,0.897076);
%                rgb(473pt)=(0.838513,0.861335,0.895991);
%                rgb(474pt)=(0.839567,0.861558,0.894901);
%                rgb(475pt)=(0.840618,0.861775,0.893805);
%                rgb(476pt)=(0.841666,0.861987,0.892704);
%                rgb(477pt)=(0.84271,0.862194,0.891599);
%                rgb(478pt)=(0.843752,0.862395,0.890488);
%                rgb(479pt)=(0.844791,0.862591,0.889372);
%                rgb(480pt)=(0.845827,0.862782,0.888251);
%                rgb(481pt)=(0.84686,0.862967,0.887125);
%                rgb(482pt)=(0.847889,0.863146,0.885994);
%                rgb(483pt)=(0.848916,0.863321,0.884858);
%                rgb(484pt)=(0.849939,0.863489,0.883717);
%                rgb(485pt)=(0.850959,0.863653,0.88257);
%                rgb(486pt)=(0.851976,0.863811,0.881419);
%                rgb(487pt)=(0.85299,0.863963,0.880263);
%                rgb(488pt)=(0.854001,0.86411,0.879102);
%                rgb(489pt)=(0.855009,0.864252,0.877936);
%                rgb(490pt)=(0.856013,0.864388,0.876765);
%                rgb(491pt)=(0.857015,0.864519,0.875589);
%                rgb(492pt)=(0.858013,0.864644,0.874408);
%                rgb(493pt)=(0.859008,0.864764,0.873222);
%                rgb(494pt)=(0.859999,0.864879,0.872032);
%                rgb(495pt)=(0.860988,0.864988,0.870836);
%                rgb(496pt)=(0.861973,0.865091,0.869636);
%                rgb(497pt)=(0.862955,0.865189,0.86843);
%                rgb(498pt)=(0.863933,0.865282,0.86722);
%                rgb(499pt)=(0.864909,0.865369,0.866005);
%                rgb(500pt)=(0.865963,0.865187,0.8647);
%                rgb(501pt)=(0.867093,0.864738,0.863308);
%                rgb(502pt)=(0.868218,0.864283,0.861914);
%                rgb(503pt)=(0.869336,0.863823,0.860517);
%                rgb(504pt)=(0.870448,0.863357,0.859118);
%                rgb(505pt)=(0.871553,0.862887,0.857717);
%                rgb(506pt)=(0.872652,0.86241,0.856313);
%                rgb(507pt)=(0.873745,0.861929,0.854907);
%                rgb(508pt)=(0.874832,0.861442,0.853499);
%                rgb(509pt)=(0.875912,0.86095,0.852088);
%                rgb(510pt)=(0.876986,0.860452,0.850676);
%                rgb(511pt)=(0.878054,0.85995,0.849261);
%                rgb(512pt)=(0.879116,0.859441,0.847843);
%                rgb(513pt)=(0.880171,0.858928,0.846424);
%                rgb(514pt)=(0.88122,0.858409,0.845002);
%                rgb(515pt)=(0.882263,0.857885,0.843579);
%                rgb(516pt)=(0.8833,0.857356,0.842153);
%                rgb(517pt)=(0.88433,0.856821,0.840725);
%                rgb(518pt)=(0.885354,0.856282,0.839295);
%                rgb(519pt)=(0.886372,0.855736,0.837862);
%                rgb(520pt)=(0.887384,0.855186,0.836428);
%                rgb(521pt)=(0.888389,0.85463,0.834992);
%                rgb(522pt)=(0.889388,0.854069,0.833553);
%                rgb(523pt)=(0.890381,0.853503,0.832113);
%                rgb(524pt)=(0.891368,0.852932,0.83067);
%                rgb(525pt)=(0.892348,0.852355,0.829225);
%                rgb(526pt)=(0.893323,0.851773,0.827779);
%                rgb(527pt)=(0.894291,0.851186,0.82633);
%                rgb(528pt)=(0.895252,0.850594,0.82488);
%                rgb(529pt)=(0.896208,0.849996,0.823427);
%                rgb(530pt)=(0.897157,0.849393,0.821973);
%                rgb(531pt)=(0.898101,0.848785,0.820517);
%                rgb(532pt)=(0.899037,0.848172,0.819059);
%                rgb(533pt)=(0.899968,0.847553,0.817598);
%                rgb(534pt)=(0.900893,0.84693,0.816136);
%                rgb(535pt)=(0.901811,0.846301,0.814673);
%                rgb(536pt)=(0.902723,0.845667,0.813207);
%                rgb(537pt)=(0.903629,0.845028,0.811739);
%                rgb(538pt)=(0.904529,0.844383,0.81027);
%                rgb(539pt)=(0.905422,0.843734,0.808799);
%                rgb(540pt)=(0.90631,0.843079,0.807326);
%                rgb(541pt)=(0.907191,0.842419,0.805851);
%                rgb(542pt)=(0.908066,0.841754,0.804375);
%                rgb(543pt)=(0.908934,0.841084,0.802897);
%                rgb(544pt)=(0.909797,0.840408,0.801417);
%                rgb(545pt)=(0.910653,0.839728,0.799935);
%                rgb(546pt)=(0.911503,0.839042,0.798452);
%                rgb(547pt)=(0.912347,0.838351,0.796967);
%                rgb(548pt)=(0.913185,0.837655,0.79548);
%                rgb(549pt)=(0.914017,0.836954,0.793992);
%                rgb(550pt)=(0.914842,0.836248,0.792502);
%                rgb(551pt)=(0.915662,0.835537,0.79101);
%                rgb(552pt)=(0.916475,0.83482,0.789517);
%                rgb(553pt)=(0.917282,0.834099,0.788022);
%                rgb(554pt)=(0.918082,0.833372,0.786526);
%                rgb(555pt)=(0.918877,0.832641,0.785028);
%                rgb(556pt)=(0.919665,0.831904,0.783528);
%                rgb(557pt)=(0.920447,0.831162,0.782027);
%                rgb(558pt)=(0.921223,0.830416,0.780525);
%                rgb(559pt)=(0.921993,0.829664,0.779021);
%                rgb(560pt)=(0.922757,0.828907,0.777515);
%                rgb(561pt)=(0.923515,0.828145,0.776008);
%                rgb(562pt)=(0.924266,0.827378,0.7745);
%                rgb(563pt)=(0.925011,0.826606,0.77299);
%                rgb(564pt)=(0.92575,0.825829,0.771479);
%                rgb(565pt)=(0.926483,0.825046,0.769966);
%                rgb(566pt)=(0.92721,0.824259,0.768452);
%                rgb(567pt)=(0.927931,0.823467,0.766937);
%                rgb(568pt)=(0.928645,0.82267,0.76542);
%                rgb(569pt)=(0.929354,0.821868,0.763902);
%                rgb(570pt)=(0.930056,0.821061,0.762382);
%                rgb(571pt)=(0.930752,0.820249,0.760861);
%                rgb(572pt)=(0.931442,0.819432,0.759339);
%                rgb(573pt)=(0.932125,0.81861,0.757815);
%                rgb(574pt)=(0.932803,0.817783,0.756291);
%                rgb(575pt)=(0.933474,0.816951,0.754765);
%                rgb(576pt)=(0.93414,0.816114,0.753237);
%                rgb(577pt)=(0.934799,0.815272,0.751709);
%                rgb(578pt)=(0.935452,0.814425,0.750179);
%                rgb(579pt)=(0.936099,0.813573,0.748648);
%                rgb(580pt)=(0.936739,0.812717,0.747116);
%                rgb(581pt)=(0.937374,0.811855,0.745583);
%                rgb(582pt)=(0.938002,0.810989,0.744048);
%                rgb(583pt)=(0.938625,0.810117,0.742513);
%                rgb(584pt)=(0.939241,0.809241,0.740976);
%                rgb(585pt)=(0.939851,0.80836,0.739438);
%                rgb(586pt)=(0.940455,0.807474,0.737899);
%                rgb(587pt)=(0.941053,0.806583,0.736359);
%                rgb(588pt)=(0.941645,0.805687,0.734818);
%                rgb(589pt)=(0.94223,0.804786,0.733276);
%                rgb(590pt)=(0.94281,0.803881,0.731732);
%                rgb(591pt)=(0.943383,0.80297,0.730188);
%                rgb(592pt)=(0.94395,0.802055,0.728643);
%                rgb(593pt)=(0.944511,0.801135,0.727097);
%                rgb(594pt)=(0.945066,0.80021,0.725549);
%                rgb(595pt)=(0.945615,0.79928,0.724001);
%                rgb(596pt)=(0.946158,0.798346,0.722452);
%                rgb(597pt)=(0.946695,0.797407,0.720902);
%                rgb(598pt)=(0.947225,0.796462,0.71935);
%                rgb(599pt)=(0.947749,0.795514,0.717798);
%                rgb(600pt)=(0.948268,0.79456,0.716245);
%                rgb(601pt)=(0.94878,0.793601,0.714692);
%                rgb(602pt)=(0.949286,0.792638,0.713137);
%                rgb(603pt)=(0.949786,0.79167,0.711581);
%                rgb(604pt)=(0.95028,0.790697,0.710025);
%                rgb(605pt)=(0.950768,0.78972,0.708468);
%                rgb(606pt)=(0.95125,0.788738,0.70691);
%                rgb(607pt)=(0.951725,0.787751,0.705351);
%                rgb(608pt)=(0.952195,0.786759,0.703791);
%                rgb(609pt)=(0.952658,0.785762,0.702231);
%                rgb(610pt)=(0.953115,0.784761,0.700669);
%                rgb(611pt)=(0.953567,0.783756,0.699107);
%                rgb(612pt)=(0.954012,0.782745,0.697545);
%                rgb(613pt)=(0.954451,0.78173,0.695981);
%                rgb(614pt)=(0.954884,0.78071,0.694417);
%                rgb(615pt)=(0.955311,0.779685,0.692852);
%                rgb(616pt)=(0.955731,0.778656,0.691287);
%                rgb(617pt)=(0.956146,0.777622,0.689721);
%                rgb(618pt)=(0.956555,0.776584,0.688154);
%                rgb(619pt)=(0.956957,0.77554,0.686587);
%                rgb(620pt)=(0.957354,0.774493,0.685019);
%                rgb(621pt)=(0.957744,0.77344,0.68345);
%                rgb(622pt)=(0.958128,0.772383,0.681881);
%                rgb(623pt)=(0.958507,0.771321,0.680311);
%                rgb(624pt)=(0.958879,0.770255,0.67874);
%                rgb(625pt)=(0.959245,0.769184,0.677169);
%                rgb(626pt)=(0.959605,0.768109,0.675598);
%                rgb(627pt)=(0.959959,0.767029,0.674026);
%                rgb(628pt)=(0.960307,0.765944,0.672453);
%                rgb(629pt)=(0.960649,0.764855,0.67088);
%                rgb(630pt)=(0.960985,0.763761,0.669307);
%                rgb(631pt)=(0.961314,0.762663,0.667733);
%                rgb(632pt)=(0.961638,0.76156,0.666158);
%                rgb(633pt)=(0.961956,0.760453,0.664583);
%                rgb(634pt)=(0.962267,0.759341,0.663008);
%                rgb(635pt)=(0.962573,0.758224,0.661432);
%                rgb(636pt)=(0.962872,0.757103,0.659856);
%                rgb(637pt)=(0.963166,0.755978,0.65828);
%                rgb(638pt)=(0.963453,0.754848,0.656703);
%                rgb(639pt)=(0.963734,0.753714,0.655125);
%                rgb(640pt)=(0.96401,0.752575,0.653548);
%                rgb(641pt)=(0.964279,0.751431,0.65197);
%                rgb(642pt)=(0.964542,0.750283,0.650391);
%                rgb(643pt)=(0.964799,0.749131,0.648813);
%                rgb(644pt)=(0.96505,0.747974,0.647234);
%                rgb(645pt)=(0.965295,0.746813,0.645655);
%                rgb(646pt)=(0.965535,0.745648,0.644075);
%                rgb(647pt)=(0.965768,0.744478,0.642496);
%                rgb(648pt)=(0.965995,0.743303,0.640916);
%                rgb(649pt)=(0.966216,0.742124,0.639335);
%                rgb(650pt)=(0.966431,0.740941,0.637755);
%                rgb(651pt)=(0.96664,0.739754,0.636174);
%                rgb(652pt)=(0.966843,0.738562,0.634594);
%                rgb(653pt)=(0.96704,0.737365,0.633013);
%                rgb(654pt)=(0.96723,0.736165,0.631431);
%                rgb(655pt)=(0.967415,0.734959,0.62985);
%                rgb(656pt)=(0.967594,0.73375,0.628269);
%                rgb(657pt)=(0.967767,0.732536,0.626687);
%                rgb(658pt)=(0.967934,0.731318,0.625105);
%                rgb(659pt)=(0.968095,0.730096,0.623524);
%                rgb(660pt)=(0.96825,0.728869,0.621942);
%                rgb(661pt)=(0.968399,0.727638,0.62036);
%                rgb(662pt)=(0.968542,0.726403,0.618778);
%                rgb(663pt)=(0.968679,0.725163,0.617196);
%                rgb(664pt)=(0.96881,0.723919,0.615614);
%                rgb(665pt)=(0.968934,0.722671,0.614032);
%                rgb(666pt)=(0.969053,0.721418,0.612449);
%                rgb(667pt)=(0.969166,0.720162,0.610867);
%                rgb(668pt)=(0.969273,0.718901,0.609285);
%                rgb(669pt)=(0.969374,0.717636,0.607703);
%                rgb(670pt)=(0.969469,0.716366,0.606121);
%                rgb(671pt)=(0.969559,0.715092,0.604539);
%                rgb(672pt)=(0.969642,0.713815,0.602957);
%                rgb(673pt)=(0.969719,0.712533,0.601375);
%                rgb(674pt)=(0.96979,0.711246,0.599794);
%                rgb(675pt)=(0.969855,0.709956,0.598212);
%                rgb(676pt)=(0.969914,0.708661,0.596631);
%                rgb(677pt)=(0.969968,0.707362,0.595049);
%                rgb(678pt)=(0.970015,0.706059,0.593468);
%                rgb(679pt)=(0.970056,0.704752,0.591887);
%                rgb(680pt)=(0.970092,0.703441,0.590306);
%                rgb(681pt)=(0.970121,0.702125,0.588725);
%                rgb(682pt)=(0.970145,0.700806,0.587145);
%                rgb(683pt)=(0.970163,0.699482,0.585564);
%                rgb(684pt)=(0.970174,0.698154,0.583984);
%                rgb(685pt)=(0.97018,0.696823,0.582404);
%                rgb(686pt)=(0.97018,0.695487,0.580825);
%                rgb(687pt)=(0.970174,0.694146,0.579245);
%                rgb(688pt)=(0.970162,0.692802,0.577666);
%                rgb(689pt)=(0.970144,0.691454,0.576087);
%                rgb(690pt)=(0.97012,0.690102,0.574509);
%                rgb(691pt)=(0.97009,0.688745,0.572931);
%                rgb(692pt)=(0.970055,0.687385,0.571353);
%                rgb(693pt)=(0.970013,0.68602,0.569775);
%                rgb(694pt)=(0.969965,0.684652,0.568198);
%                rgb(695pt)=(0.969912,0.683279,0.566621);
%                rgb(696pt)=(0.969853,0.681903,0.565044);
%                rgb(697pt)=(0.969787,0.680522,0.563468);
%                rgb(698pt)=(0.969716,0.679138,0.561892);
%                rgb(699pt)=(0.969639,0.677749,0.560317);
%                rgb(700pt)=(0.969557,0.676357,0.558742);
%                rgb(701pt)=(0.969468,0.67496,0.557167);
%                rgb(702pt)=(0.969373,0.67356,0.555593);
%                rgb(703pt)=(0.969273,0.672155,0.55402);
%                rgb(704pt)=(0.969166,0.670747,0.552446);
%                rgb(705pt)=(0.969054,0.669335,0.550874);
%                rgb(706pt)=(0.968936,0.667918,0.549301);
%                rgb(707pt)=(0.968812,0.666498,0.547729);
%                rgb(708pt)=(0.968682,0.665074,0.546158);
%                rgb(709pt)=(0.968546,0.663646,0.544587);
%                rgb(710pt)=(0.968405,0.662214,0.543017);
%                rgb(711pt)=(0.968257,0.660778,0.541448);
%                rgb(712pt)=(0.968104,0.659339,0.539878);
%                rgb(713pt)=(0.967945,0.657895,0.53831);
%                rgb(714pt)=(0.96778,0.656448,0.536742);
%                rgb(715pt)=(0.967609,0.654996,0.535174);
%                rgb(716pt)=(0.967433,0.653541,0.533608);
%                rgb(717pt)=(0.96725,0.652082,0.532041);
%                rgb(718pt)=(0.967062,0.650619,0.530476);
%                rgb(719pt)=(0.966868,0.649153,0.528911);
%                rgb(720pt)=(0.966668,0.647682,0.527347);
%                rgb(721pt)=(0.966463,0.646208,0.525783);
%                rgb(722pt)=(0.966251,0.64473,0.52422);
%                rgb(723pt)=(0.966034,0.643248,0.522658);
%                rgb(724pt)=(0.965811,0.641762,0.521096);
%                rgb(725pt)=(0.965582,0.640273,0.519535);
%                rgb(726pt)=(0.965347,0.638779,0.517975);
%                rgb(727pt)=(0.965107,0.637282,0.516415);
%                rgb(728pt)=(0.964861,0.635782,0.514857);
%                rgb(729pt)=(0.964609,0.634277,0.513299);
%                rgb(730pt)=(0.964351,0.632769,0.511741);
%                rgb(731pt)=(0.964087,0.631257,0.510185);
%                rgb(732pt)=(0.963818,0.629741,0.508629);
%                rgb(733pt)=(0.963543,0.628222,0.507074);
%                rgb(734pt)=(0.963262,0.626698,0.50552);
%                rgb(735pt)=(0.962975,0.625171,0.503967);
%                rgb(736pt)=(0.962683,0.623641,0.502414);
%                rgb(737pt)=(0.962385,0.622106,0.500863);
%                rgb(738pt)=(0.962081,0.620569,0.499312);
%                rgb(739pt)=(0.961772,0.619027,0.497762);
%                rgb(740pt)=(0.961456,0.617482,0.496213);
%                rgb(741pt)=(0.961135,0.615933,0.494664);
%                rgb(742pt)=(0.960809,0.61438,0.493117);
%                rgb(743pt)=(0.960476,0.612824,0.491571);
%                rgb(744pt)=(0.960138,0.611264,0.490025);
%                rgb(745pt)=(0.959794,0.6097,0.488481);
%                rgb(746pt)=(0.959444,0.608133,0.486937);
%                rgb(747pt)=(0.959089,0.606562,0.485394);
%                rgb(748pt)=(0.958728,0.604987,0.483852);
%                rgb(749pt)=(0.958362,0.603409,0.482311);
%                rgb(750pt)=(0.957989,0.601828,0.480772);
%                rgb(751pt)=(0.957611,0.600242,0.479233);
%                rgb(752pt)=(0.957227,0.598654,0.477695);
%                rgb(753pt)=(0.956838,0.597061,0.476158);
%                rgb(754pt)=(0.956443,0.595465,0.474622);
%                rgb(755pt)=(0.956042,0.593866,0.473087);
%                rgb(756pt)=(0.955636,0.592262,0.471553);
%                rgb(757pt)=(0.955224,0.590656,0.470021);
%                rgb(758pt)=(0.954806,0.589045,0.468489);
%                rgb(759pt)=(0.954383,0.587432,0.466958);
%                rgb(760pt)=(0.953954,0.585814,0.465429);
%                rgb(761pt)=(0.95352,0.584193,0.4639);
%                rgb(762pt)=(0.953079,0.582569,0.462373);
%                rgb(763pt)=(0.952634,0.580941,0.460846);
%                rgb(764pt)=(0.952182,0.57931,0.459321);
%                rgb(765pt)=(0.951725,0.577675,0.457797);
%                rgb(766pt)=(0.951262,0.576036,0.456274);
%                rgb(767pt)=(0.950794,0.574394,0.454752);
%                rgb(768pt)=(0.95032,0.572749,0.453232);
%                rgb(769pt)=(0.949841,0.5711,0.451712);
%                rgb(770pt)=(0.949356,0.569447,0.450194);
%                rgb(771pt)=(0.948865,0.567791,0.448677);
%                rgb(772pt)=(0.948369,0.566132,0.447161);
%                rgb(773pt)=(0.947867,0.564469,0.445646);
%                rgb(774pt)=(0.94736,0.562803,0.444133);
%                rgb(775pt)=(0.946847,0.561133,0.44262);
%                rgb(776pt)=(0.946329,0.55946,0.441109);
%                rgb(777pt)=(0.945805,0.557783,0.439599);
%                rgb(778pt)=(0.945276,0.556103,0.438091);
%                rgb(779pt)=(0.944741,0.554419,0.436583);
%                rgb(780pt)=(0.9442,0.552732,0.435077);
%                rgb(781pt)=(0.943654,0.551042,0.433573);
%                rgb(782pt)=(0.943102,0.549348,0.432069);
%                rgb(783pt)=(0.942545,0.547651,0.430567);
%                rgb(784pt)=(0.941983,0.54595,0.429066);
%                rgb(785pt)=(0.941415,0.544246,0.427567);
%                rgb(786pt)=(0.940841,0.542538,0.426068);
%                rgb(787pt)=(0.940262,0.540827,0.424572);
%                rgb(788pt)=(0.939677,0.539113,0.423076);
%                rgb(789pt)=(0.939087,0.537395,0.421582);
%                rgb(790pt)=(0.938492,0.535674,0.420089);
%                rgb(791pt)=(0.937891,0.533949,0.418598);
%                rgb(792pt)=(0.937285,0.532221,0.417108);
%                rgb(793pt)=(0.936673,0.53049,0.415619);
%                rgb(794pt)=(0.936055,0.528755,0.414132);
%                rgb(795pt)=(0.935433,0.527017,0.412646);
%                rgb(796pt)=(0.934804,0.525276,0.411161);
%                rgb(797pt)=(0.934171,0.523531,0.409678);
%                rgb(798pt)=(0.933532,0.521782,0.408197);
%                rgb(799pt)=(0.932887,0.520031,0.406717);
%                rgb(800pt)=(0.932238,0.518276,0.405238);
%                rgb(801pt)=(0.931582,0.516517,0.403761);
%                rgb(802pt)=(0.930922,0.514755,0.402285);
%                rgb(803pt)=(0.930256,0.51299,0.400811);
%                rgb(804pt)=(0.929584,0.511222,0.399338);
%                rgb(805pt)=(0.928907,0.50945,0.397867);
%                rgb(806pt)=(0.928225,0.507674,0.396397);
%                rgb(807pt)=(0.927538,0.505896,0.394929);
%                rgb(808pt)=(0.926845,0.504114,0.393462);
%                rgb(809pt)=(0.926147,0.502328,0.391997);
%                rgb(810pt)=(0.925443,0.500539,0.390533);
%                rgb(811pt)=(0.924734,0.498747,0.389071);
%                rgb(812pt)=(0.92402,0.496952,0.387611);
%                rgb(813pt)=(0.9233,0.495153,0.386152);
%                rgb(814pt)=(0.922576,0.49335,0.384694);
%                rgb(815pt)=(0.921845,0.491545,0.383239);
%                rgb(816pt)=(0.92111,0.489736,0.381784);
%                rgb(817pt)=(0.920369,0.487923,0.380332);
%                rgb(818pt)=(0.919623,0.486107,0.378881);
%                rgb(819pt)=(0.918872,0.484288,0.377432);
%                rgb(820pt)=(0.918115,0.482465,0.375984);
%                rgb(821pt)=(0.917353,0.480639,0.374538);
%                rgb(822pt)=(0.916586,0.47881,0.373093);
%                rgb(823pt)=(0.915814,0.476977,0.37165);
%                rgb(824pt)=(0.915036,0.475141,0.370209);
%                rgb(825pt)=(0.914253,0.473301,0.36877);
%                rgb(826pt)=(0.913465,0.471458,0.367332);
%                rgb(827pt)=(0.912672,0.469612,0.365896);
%                rgb(828pt)=(0.911874,0.467762,0.364462);
%                rgb(829pt)=(0.91107,0.465909,0.363029);
%                rgb(830pt)=(0.910261,0.464052,0.361598);
%                rgb(831pt)=(0.909447,0.462192,0.360169);
%                rgb(832pt)=(0.908627,0.460328,0.358741);
%                rgb(833pt)=(0.907803,0.458461,0.357315);
%                rgb(834pt)=(0.906973,0.45659,0.355891);
%                rgb(835pt)=(0.906138,0.454716,0.354469);
%                rgb(836pt)=(0.905298,0.452839,0.353048);
%                rgb(837pt)=(0.904453,0.450958,0.35163);
%                rgb(838pt)=(0.903603,0.449073,0.350213);
%                rgb(839pt)=(0.902748,0.447185,0.348797);
%                rgb(840pt)=(0.901887,0.445294,0.347384);
%                rgb(841pt)=(0.901021,0.443399,0.345972);
%                rgb(842pt)=(0.900151,0.4415,0.344562);
%                rgb(843pt)=(0.899275,0.439598,0.343154);
%                rgb(844pt)=(0.898394,0.437693,0.341748);
%                rgb(845pt)=(0.897508,0.435784,0.340344);
%                rgb(846pt)=(0.896617,0.433871,0.338941);
%                rgb(847pt)=(0.89572,0.431954,0.33754);
%                rgb(848pt)=(0.894819,0.430035,0.336141);
%                rgb(849pt)=(0.893913,0.428111,0.334744);
%                rgb(850pt)=(0.893002,0.426184,0.333349);
%                rgb(851pt)=(0.892085,0.424253,0.331956);
%                rgb(852pt)=(0.891164,0.422318,0.330565);
%                rgb(853pt)=(0.890237,0.42038,0.329175);
%                rgb(854pt)=(0.889306,0.418438,0.327788);
%                rgb(855pt)=(0.888369,0.416493,0.326402);
%                rgb(856pt)=(0.887428,0.414543,0.325018);
%                rgb(857pt)=(0.886481,0.41259,0.323636);
%                rgb(858pt)=(0.88553,0.410634,0.322256);
%                rgb(859pt)=(0.884574,0.408673,0.320878);
%                rgb(860pt)=(0.883612,0.406709,0.319502);
%                rgb(861pt)=(0.882646,0.40474,0.318128);
%                rgb(862pt)=(0.881674,0.402768,0.316756);
%                rgb(863pt)=(0.880698,0.400792,0.315385);
%                rgb(864pt)=(0.879717,0.398812,0.314017);
%                rgb(865pt)=(0.878731,0.396829,0.312651);
%                rgb(866pt)=(0.87774,0.394841,0.311287);
%                rgb(867pt)=(0.876744,0.392849,0.309924);
%                rgb(868pt)=(0.875743,0.390854,0.308564);
%                rgb(869pt)=(0.874738,0.388854,0.307206);
%                rgb(870pt)=(0.873727,0.38685,0.305849);
%                rgb(871pt)=(0.872712,0.384842,0.304495);
%                rgb(872pt)=(0.871691,0.38283,0.303143);
%                rgb(873pt)=(0.870666,0.380814,0.301792);
%                rgb(874pt)=(0.869636,0.378794,0.300444);
%                rgb(875pt)=(0.868601,0.37677,0.299098);
%                rgb(876pt)=(0.867562,0.374741,0.297754);
%                rgb(877pt)=(0.866517,0.372708,0.296412);
%                rgb(878pt)=(0.865468,0.370671,0.295072);
%                rgb(879pt)=(0.864414,0.368629,0.293734);
%                rgb(880pt)=(0.863355,0.366583,0.292398);
%                rgb(881pt)=(0.862291,0.364533,0.291064);
%                rgb(882pt)=(0.861223,0.362478,0.289733);
%                rgb(883pt)=(0.860149,0.360418,0.288403);
%                rgb(884pt)=(0.859071,0.358354,0.287075);
%                rgb(885pt)=(0.857989,0.356286,0.28575);
%                rgb(886pt)=(0.856901,0.354213,0.284427);
%                rgb(887pt)=(0.855809,0.352135,0.283106);
%                rgb(888pt)=(0.854712,0.350052,0.281787);
%                rgb(889pt)=(0.853611,0.347964,0.28047);
%                rgb(890pt)=(0.852504,0.345872,0.279155);
%                rgb(891pt)=(0.851393,0.343775,0.277842);
%                rgb(892pt)=(0.850278,0.341673,0.276532);
%                rgb(893pt)=(0.849157,0.339565,0.275224);
%                rgb(894pt)=(0.848032,0.337453,0.273918);
%                rgb(895pt)=(0.846902,0.335335,0.272614);
%                rgb(896pt)=(0.845768,0.333213,0.271312);
%                rgb(897pt)=(0.844629,0.331085,0.270012);
%                rgb(898pt)=(0.843486,0.328951,0.268715);
%                rgb(899pt)=(0.842337,0.326812,0.26742);
%                rgb(900pt)=(0.841185,0.324668,0.266127);
%                rgb(901pt)=(0.840027,0.322518,0.264836);
%                rgb(902pt)=(0.838865,0.320363,0.263548);
%                rgb(903pt)=(0.837699,0.318201,0.262261);
%                rgb(904pt)=(0.836528,0.316034,0.260977);
%                rgb(905pt)=(0.835352,0.313861,0.259695);
%                rgb(906pt)=(0.834172,0.311682,0.258416);
%                rgb(907pt)=(0.832987,0.309497,0.257138);
%                rgb(908pt)=(0.831798,0.307306,0.255863);
%                rgb(909pt)=(0.830604,0.305108,0.25459);
%                rgb(910pt)=(0.829406,0.302904,0.25332);
%                rgb(911pt)=(0.828203,0.300694,0.252052);
%                rgb(912pt)=(0.826995,0.298476,0.250785);
%                rgb(913pt)=(0.825784,0.296253,0.249522);
%                rgb(914pt)=(0.824567,0.294022,0.24826);
%                rgb(915pt)=(0.823347,0.291784,0.247001);
%                rgb(916pt)=(0.822122,0.289539,0.245744);
%                rgb(917pt)=(0.820892,0.287287,0.24449);
%                rgb(918pt)=(0.819658,0.285028,0.243237);
%                rgb(919pt)=(0.81842,0.282761,0.241988);
%                rgb(920pt)=(0.817177,0.280486,0.24074);
%                rgb(921pt)=(0.815929,0.278204,0.239495);
%                rgb(922pt)=(0.814678,0.275913,0.238252);
%                rgb(923pt)=(0.813422,0.273615,0.237011);
%                rgb(924pt)=(0.812161,0.271308,0.235773);
%                rgb(925pt)=(0.810897,0.268993,0.234537);
%                rgb(926pt)=(0.809628,0.266669,0.233304);
%                rgb(927pt)=(0.808354,0.264336,0.232072);
%                rgb(928pt)=(0.807077,0.261994,0.230844);
%                rgb(929pt)=(0.805795,0.259642,0.229617);
%                rgb(930pt)=(0.804508,0.257282,0.228393);
%                rgb(931pt)=(0.803218,0.254911,0.227171);
%                rgb(932pt)=(0.801923,0.252531,0.225952);
%                rgb(933pt)=(0.800624,0.25014,0.224735);
%                rgb(934pt)=(0.79932,0.247739,0.223521);
%                rgb(935pt)=(0.798013,0.245327,0.222309);
%                rgb(936pt)=(0.796701,0.242904,0.221099);
%                rgb(937pt)=(0.795384,0.24047,0.219892);
%                rgb(938pt)=(0.794064,0.238024,0.218687);
%                rgb(939pt)=(0.79274,0.235566,0.217485);
%                rgb(940pt)=(0.791411,0.233096,0.216285);
%                rgb(941pt)=(0.790078,0.230614,0.215088);
%                rgb(942pt)=(0.788741,0.228118,0.213893);
%                rgb(943pt)=(0.787399,0.225609,0.2127);
%                rgb(944pt)=(0.786054,0.223086,0.21151);
%                rgb(945pt)=(0.784705,0.220549,0.210322);
%                rgb(946pt)=(0.783351,0.217998,0.209137);
%                rgb(947pt)=(0.781993,0.215431,0.207955);
%                rgb(948pt)=(0.780631,0.212849,0.206774);
%                rgb(949pt)=(0.779265,0.21025,0.205597);
%                rgb(950pt)=(0.777895,0.207635,0.204421);
%                rgb(951pt)=(0.776521,0.205003,0.203249);
%                rgb(952pt)=(0.775143,0.202353,0.202078);
%                rgb(953pt)=(0.77376,0.199685,0.200911);
%                rgb(954pt)=(0.772374,0.196997,0.199746);
%                rgb(955pt)=(0.770984,0.194289,0.198583);
%                rgb(956pt)=(0.769589,0.191561,0.197423);
%                rgb(957pt)=(0.768191,0.188812,0.196265);
%                rgb(958pt)=(0.766789,0.18604,0.19511);
%                rgb(959pt)=(0.765382,0.183245,0.193957);
%                rgb(960pt)=(0.763972,0.180426,0.192807);
%                rgb(961pt)=(0.762558,0.177582,0.19166);
%                rgb(962pt)=(0.761139,0.174711,0.190515);
%                rgb(963pt)=(0.759717,0.171813,0.189373);
%                rgb(964pt)=(0.758291,0.168887,0.188233);
%                rgb(965pt)=(0.756861,0.16593,0.187096);
%                rgb(966pt)=(0.755427,0.162942,0.185961);
%                rgb(967pt)=(0.753989,0.15992,0.184829);
%                rgb(968pt)=(0.752548,0.156864,0.1837);
%                rgb(969pt)=(0.751102,0.153771,0.182573);
%                rgb(970pt)=(0.749653,0.15064,0.181449);
%                rgb(971pt)=(0.7482,0.147468,0.180327);
%                rgb(972pt)=(0.746743,0.144253,0.179208);
%                rgb(973pt)=(0.745282,0.140992,0.178092);
%                rgb(974pt)=(0.743817,0.137682,0.176978);
%                rgb(975pt)=(0.742349,0.134321,0.175867);
%                rgb(976pt)=(0.740876,0.130906,0.174758);
%                rgb(977pt)=(0.7394,0.127431,0.173652);
%                rgb(978pt)=(0.73792,0.123894,0.172549);
%                rgb(979pt)=(0.736437,0.120289,0.171448);
%                rgb(980pt)=(0.73495,0.116611,0.170351);
%                rgb(981pt)=(0.733459,0.112854,0.169255);
%                rgb(982pt)=(0.731964,0.109013,0.168163);
%                rgb(983pt)=(0.730465,0.105078,0.167073);
%                rgb(984pt)=(0.728963,0.101042,0.165985);
%                rgb(985pt)=(0.727458,0.0968943,0.164901);
%                rgb(986pt)=(0.725948,0.0926233,0.163819);
%                rgb(987pt)=(0.724435,0.0882153,0.16274);
%                rgb(988pt)=(0.722918,0.083654,0.161663);
%                rgb(989pt)=(0.721398,0.0789198,0.160589);
%                rgb(990pt)=(0.719874,0.0739888,0.159518);
%                rgb(991pt)=(0.718346,0.0688314,0.15845);
%                rgb(992pt)=(0.716815,0.0634104,0.157384);
%                rgb(993pt)=(0.715281,0.0576774,0.156321);
%                rgb(994pt)=(0.713742,0.0515686,0.155261);
%                rgb(995pt)=(0.7122,0.0449959,0.154204);
%                rgb(996pt)=(0.710655,0.037839,0.153149);
%                rgb(997pt)=(0.709106,0.0305255,0.152097);
%                rgb(998pt)=(0.707554,0.0232881,0.151048);
%                rgb(999pt)=(0.705998,0.0161265,0.150001)}
%%

DGforM1.tex
\documentclass[final,preprint,3p,sort&compress]{elsarticle}
%% to remove the Elsevier thing on the first page
\journal{{\tt arXiv.org}}
%\journal{{\tt Journal of Computational and Applied Mathematics}}
\usepackage[dvips]{epsfig}
\usepackage{graphicx} 
\usepackage{pdfpages}
\usepackage{latexsym}
\usepackage{verbatim}
\usepackage{doi}        
\usepackage{amsmath,amsthm}
%\usepackage{bm}
%\usepackage{amsthm}
\usepackage{amssymb}
%\usepackage{a4wide}
\usepackage{bbm}
%\usepackage{mathtools}
%\usepackage{pgfplots}
%\pgfplotsset{compat=newest}
%\usepackage{refcheck}
%\usepackage{fonts}
\usepackage{enumerate}
%\usepackage[sort]{cite}
\usepackage{placeins}
\usepackage{standalone}
\usepackage{paralist}
\usepackage{subcaption}
%\usepackage[top=1in, bottom=1.25in, left=1.0in, right=1.0in]{geometry}
\usepackage{booktabs} % Required for nicer horizontal rules in tables
\usepackage{standalone}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\definecolor{greenyellow}   {cmyk}{0.15, 0   , 0.69, 0   }
\definecolor{yellow}        {cmyk}{0   , 0   , 1   , 0   }
\definecolor{goldenrod}     {cmyk}{0   , 0.10, 0.84, 0   }
\definecolor{dandelion}     {cmyk}{0   , 0.29, 0.84, 0   }
\definecolor{apricot}       {cmyk}{0   , 0.32, 0.52, 0   }
\definecolor{peach}         {cmyk}{0   , 0.50, 0.70, 0   }
\definecolor{melon}         {cmyk}{0   , 0.46, 0.50, 0   }
\definecolor{yelloworange}  {cmyk}{0   , 0.42, 1   , 0   }
\definecolor{orange}        {cmyk}{0   , 0.61, 0.87, 0   }
\definecolor{burntorange}   {cmyk}{0   , 0.51, 1   , 0   }
\definecolor{bittersweet}   {cmyk}{0   , 0.75, 1   , 0.24}
\definecolor{redorange}     {cmyk}{0   , 0.77, 0.87, 0   }
\definecolor{mahogany}      {cmyk}{0   , 0.85, 0.87, 0.35}
\definecolor{maroon}        {cmyk}{0   , 0.87, 0.68, 0.32}
\definecolor{brickred}      {cmyk}{0   , 0.89, 0.94, 0.28}
\definecolor{red}           {cmyk}{0   , 1   , 1   , 0   }
\definecolor{orangered}     {cmyk}{0   , 1   , 0.50, 0   }
\definecolor{rubinered}     {cmyk}{0   , 1   , 0.13, 0   }
\definecolor{wildstrawberry}{cmyk}{0   , 0.96, 0.39, 0   }
\definecolor{salmon}        {cmyk}{0   , 0.53, 0.38, 0   }
\definecolor{carnationpink} {cmyk}{0   , 0.63, 0   , 0   }
\definecolor{magenta}       {cmyk}{0   , 1   , 0   , 0   }
\definecolor{violetred}     {cmyk}{0   , 0.81, 0   , 0   }
\definecolor{rhodamine}     {cmyk}{0   , 0.82, 0   , 0   }
\definecolor{mulberry}      {cmyk}{0.34, 0.90, 0   , 0.02}
\definecolor{redviolet}     {cmyk}{0.07, 0.90, 0   , 0.34}
\definecolor{fuchsia}       {cmyk}{0.47, 0.91, 0   , 0.08}
\definecolor{lavender}      {cmyk}{0   , 0.48, 0   , 0   }
\definecolor{thistle}       {cmyk}{0.12, 0.59, 0   , 0   }
\definecolor{orchid}        {cmyk}{0.32, 0.64, 0   , 0   }
\definecolor{darkorchid}    {cmyk}{0.40, 0.80, 0.20, 0   }
\definecolor{purple}        {cmyk}{0.45, 0.86, 0   , 0   }
\definecolor{plum}          {cmyk}{0.50, 1   , 0   , 0   }
\definecolor{violet}        {cmyk}{0.79, 0.88, 0   , 0   }
\definecolor{royalpurple}   {cmyk}{0.75, 0.90, 0   , 0   }
\definecolor{blueviolet}    {cmyk}{0.86, 0.91, 0   , 0.04}
\definecolor{periwinkle}    {cmyk}{0.57, 0.55, 0   , 0   }
\definecolor{cadetblue}     {cmyk}{0.62, 0.57, 0.23, 0   }
\definecolor{cornflowerblue}{cmyk}{0.65, 0.13, 0   , 0   }
\definecolor{midnightblue}  {cmyk}{0.98, 0.13, 0   , 0.43}
\definecolor{navyblue}      {cmyk}{0.94, 0.54, 0   , 0   }
\definecolor{royalblue}     {cmyk}{1   , 0.50, 0   , 0   }
\definecolor{blue}          {cmyk}{1   , 1   , 0   , 0   }
\definecolor{cerulean}      {cmyk}{0.94, 0.11, 0   , 0   }
\definecolor{cyan}          {cmyk}{1   , 0   , 0   , 0   }
\definecolor{processblue}   {cmyk}{0.96, 0   , 0   , 0   }
\definecolor{skyblue}       {cmyk}{0.62, 0   , 0.12, 0   }
\definecolor{turquoise}     {cmyk}{0.85, 0   , 0.20, 0   }
\definecolor{tealblue}      {cmyk}{0.86, 0   , 0.34, 0.02}
\definecolor{aquamarine}    {cmyk}{0.82, 0   , 0.30, 0   }
\definecolor{bluegreen}     {cmyk}{0.85, 0   , 0.33, 0   }
\definecolor{emerald}       {cmyk}{1   , 0   , 0.50, 0   }
\definecolor{junglegreen}   {cmyk}{0.99, 0   , 0.52, 0   }
\definecolor{seagreen}      {cmyk}{0.69, 0   , 0.50, 0   }
\definecolor{green}         {cmyk}{1   , 0   , 1   , 0   }
\definecolor{forestgreen}   {cmyk}{0.91, 0   , 0.88, 0.12}
\definecolor{pinegreen}     {cmyk}{0.92, 0   , 0.59, 0.25}
\definecolor{limegreen}     {cmyk}{0.50, 0   , 1   , 0   }
\definecolor{yellowgreen}   {cmyk}{0.44, 0   , 0.74, 0   }
\definecolor{springgreen}   {cmyk}{0.26, 0   , 0.76, 0   }
\definecolor{olivegreen}    {cmyk}{0.64, 0   , 0.95, 0.40}
\definecolor{rawsienna}     {cmyk}{0   , 0.72, 1   , 0.45}
\definecolor{sepia}         {cmyk}{0   , 0.83, 1   , 0.70}
\definecolor{brown}         {cmyk}{0   , 0.81, 1   , 0.60}
\definecolor{tan}           {cmyk}{0.14, 0.42, 0.56, 0   }
\definecolor{gray}          {cmyk}{0   , 0   , 0   , 0.50}
\definecolor{black}         {cmyk}{0   , 0   , 0   , 1   }
\definecolor{white}         {cmyk}{0   , 0   , 0   , 0   } 

% *************** Enable hyperlinks in PDF documents ***************
 \usepackage[]{hyperref} 

%\newif\ifPrint
%\Printtrue
%\Printfalse
        
%\ifPrint
\newcommand{\externaltikz}[2]{\includegraphics{Externals/t#1}}
%\else      
%\usepackage[ComplainAboutUnexternalized]{externaltikz}
%\fi
\usepackage{tikz,pgfplots}
\pgfplotsset{compat=1.12}
\usepgfplotslibrary{groupplots}
\input{colormap}
\newlength\figurewidth
\setlength{\figurewidth}{3.5cm}
\newlength\figureheight
\setlength{\figureheight}{\figurewidth}

\usepackage{relinput}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}[theorem]{Example}
\newtheorem{assumption}[theorem]{Assumption}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}


% ******************************************************************************
% ************************* Tikz subfigure reference ******************************
% ******************************************************************************
\newcounter{tikzsubfigcounter}[figure]
\renewcommand{\thetikzsubfigcounter}{\thesection.\the\numexpr\value{figure}+1\relax\alph{tikzsubfigcounter}}

\newcommand{\tikztitle}[1]{ %
\refstepcounter{tikzsubfigcounter}
\textbf{(\alph{tikzsubfigcounter})}\space\space #1 
}

\newcounter{tikzsubfigcounterinvisible}[figure]
\renewcommand{\thetikzsubfigcounterinvisible}{\thesection.\the\numexpr\value{figure}+1\relax\alph{tikzsubfigcounterinvisible}}

\newcommand{\settikzlabel}[1]{ %
\refstepcounter{tikzsubfigcounterinvisible} \label{#1} 
}



\newtheorem{thm}{\bf Theorem}[section]
\newtheorem{lem}[thm]{\bf Lemma}
\newtheorem{prop}[thm]{\bf Proposition}
\newtheorem{algo}[thm]{\bf Algorithm}
\newtheorem{ass}{\bf Assumption}
\newtheorem{cor}[thm]{\bf Corollary}
%\newtheorem{example}[thm]{\bf Example}

%\theoremstyle{remark}

\newtheorem{rem}[thm]{\bf  Remark}


% Befehle fuer Referenzen

\numberwithin{equation}{section}
%\numberwithin{figure}{section}
%\numberwithin{subsection}{section}
%\usepackage{fonts}
%\input{usermacros}

%--------------------------------------------------------------------------
\newcommand{\bdm}{\begin{displaymath}}
\newcommand{\edm}{\end{displaymath}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\beqa}{\begin{eqnarray}}
\newcommand{\eeqa}{\end{eqnarray}}
\parskip1.5ex
\parindent0ex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{A Comparative Study of Limiting Strategies in Discontinuous Galerkin Schemes for the $M_1$ Model of Radiation Transport}
\author[pc]{Prince Chidyagwai}
\author[mf]{Martin Frank}
\author[fs]{Florian Schneider}
\author[bs]{Benjamin Seibold}
\address[pc]{Department of Mathematics and Statistics, Loyola University Maryland, 4501 N Charles Street, Baltimore, MD 21210 {\tt pchidyagwai@loyola.edu}}
\address[mf]{Department of Mathematics, RWTH Aachen University, Schinkelstr. 2,52062 Aachen, Germany, {\tt frank@mathcces.rwth-aachen.de}}
\address[fs]{Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schr\"odinger-Str., 67663 Kaiserslautern, Germany, {\tt schneider@mathematik.uni-kl.de}}
\address[bs]{Department of Mathematics, Temple University, 1805 N Broad Street, Philadelphia, PA 19122, {\tt seibold@temple.edu}}
\date{}

\input{usermacros}

\begin{document}

\begin{abstract}
    The $M_1$ minimum entropy moment system is a system of hyperbolic balance laws that approximates the radiation transport equation, 
    and has many desirable properties. 
    Among them are symmetric hyperbolicity, entropy decay, moment realizability, and correct behavior in
    the diffusion and free-streaming limits. However, numerical difficulties arise when approximating the solution
    of the $M_1$ model by high order numerical schemes; namely maintaining the realizability of the numerical solution
    and controlling spurious oscillations. In this paper, we extend a previously constructed one-dimensional realizability limiting strategy to 2D. 
    In addition, we perform a numerical study of various combinations of the realizability limiter and the TVBM local slope limiter 
    on a third order Discontinuous Galerkin (DG) scheme on both triangular and rectangular meshes. In several test cases, 
    we demonstrate that in general, a combination of the realizability limiter and a TVBM limiter is necessary to obtain a robust and accurate numerical scheme.
    Our code is published so that all results can be reproduced by the reader.
\end{abstract}
\begin{keyword}
discontinuous Galerkin \sep moment models \sep minimum entropy \sep realizability limiter
\MSC[2010] 35L40 \sep 65M08 \sep 65M60 \sep 65M70 
\end{keyword}
\maketitle



\noindent

\input{Sections/introduction}
\input{Sections/properties}
\input{Sections/dg}
\input{Sections/realizabilitylimiter}
\input{Sections/numerics}
\input{Sections/outlook}

\section*{Acknowledgements}
%
B.\ Seibold wishes to acknowledge support by the National Science
Foundation under grant DMS-1719640. M.\ Frank acknowledges support from Deutsche Forschungsgemeinschaft (DFG) under Grant FR 2841/6-1.
We thank Philipp Monreal for co-developing the unstructured-mesh DG code.

% Bibliography
%%%%%%%%%%%%%%
\bibliographystyle{bibstyle}
\bibliography{RadLit}

\end{document}



additionalnumerics.tex
%\section{Implementation Details}
%\label{sec:details}
%
\section{Numerical Results}
\label{sec:numerical_results}
In this section, we compare the performance of the RKDG schemes on unstructured triangular meshes
and uniform rectangular meshes for various combinations of limiters. Specifically we consider the limiter combinations given in Table~\ref{tab:LimiterCombinations}.
\begin{table}[h]
\centering
\begin{tabular}{l l}
Abbr. & Explanation\\
\midrule
SL$M$$\square/\triangle$ & Slope limiter in primitive variables with the constant $M$ specified\\
CL$M$$\square/\triangle$ & Slope limiter in characteristic variables with the constant $M$ specified\\
SRL$M$$\square/\triangle$ & Slope limiter in primitive variables + realizability limiter\\
CRL$M$$\square/\triangle$ & Slope limiter in characteristic variables + realizability limiter\\
\end{tabular}
\caption{Abbreviation of limiter combinations. Setting $M=\infty$ is equivalent to disabling the slope limiter. $\square$ and $\triangle$ correspond to the rectangular and triangular meshes, respectively.}
\label{tab:LimiterCombinations}
\end{table}

A complete characterization and discussion of the first 4 test cases is presented
in the paper. The last test case is an additional test case on the collision of two beams.

\subsection{Line Source}
\begin{itemize}
\item Domain: $[-0.5,0.5]^2$ with $h =0.004$
\item Time $T=0.45$
\item Parameters: $\sigma_a=\sigma_s=0$, $q^{(0)}=0$, $q^{(1)}=0$
\item Initial condition: $\psi^{(0)}(x,y) = \max(\exp(-10\frac{x^2+y^2}{\sigma^2}),10^{-4})$ with $\sigma=0.02$, $\psi^{(1)}(x,y)=0$
\item Boundary conditions are Dirichlet conditions equal to the initial conditions. Because the signal does not reach the boundary, these conditions do not influence the result.
\end{itemize}


\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{LinesourceAdd3}{\input{Images/LinesourceAdd3}}
\caption{Line source: comparison of unlimited solution with and without realizability limiter.}
\label{fig:linesource3}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{LinesourceAdd}{\input{Images/LinesourceAdd1}}
\caption{Line source: comparison of primitive and characteristic slope limiter.}
\label{fig:linesource}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{LinesourceAdd2}{\input{Images/LinesourceAdd2}}
\caption{Line source: comparison of primitive and characteristic slope limiter with realizability limiter.}
\label{fig:linesource2}
\end{figure}


\clearpage
\subsection{Homogeneous Disk}
\begin{itemize}
\item Domain: $[-5,5]^2$ with $h=0.05$
\item Time: 3.75
\item Parameters: Let $D$ be the unit disk $D=\{ (x,y):\ x^2+y^2\leq 1\}$. Set $\sigma_s=0$; $\sigma_a=10$ on $D$, zero otherwise; source $q^{(0)}=1$ on $D$, zero otherwise; $q^{(1)} = 0$
\item Initial condition: $\psi^{(0)}=0$
\item Boundary conditions: vacuum
\end{itemize}

\begin{figure}[h]
\def\localpath{./Images/HomDisk/}
\centering
\externaltikz{HomDiskAdd3}{\input{Images/HomDiskAdd3}}
\caption{Homogeneous disk: comparison of unlimited solution with and without realizability limiter.}
\label{fig:homdisk3}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/HomDisk/}
\centering
\externaltikz{HomDiskAdd}{\input{Images/HomDiskAdd1}}
\caption{Homogeneous disk: comparison of primitive and characteristic slope limiter.}
\label{fig:homdisk}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/HomDisk/}
\centering
\externaltikz{HomDiskAdd2}{\input{Images/HomDiskAdd2}}
\caption{Homogeneous disk: comparison of primitive and characteristic slope limiter with realizability limiter.}
\label{fig:homdisk2}
\end{figure}



\clearpage
\subsection{Flash}
\begin{itemize}
\item Domain: $[-10,10]^2$ with $h =0.06$
\item Time: 6
\item Parameters: Let $D$ be unit disk with radius five $D=\{ (x,y):\ x^2+y^2\leq 25\}$. Set $\sigma_s=0$; $\sigma_a=0$; source $q^{(0)}=0$, $q^{(1)} = 0$
\item Initial condition: $\psi^{(0)}=1$, $\psi^{(1)}_x=0.9$, $\psi^{(1)}_y=0$ on $D$, zero otherwise
\item Boundary conditions: vacuum
\end{itemize}

\begin{figure}[h]
\def\localpath{./Images/Flash/}
\centering
\externaltikz{FlashAdd3}{\input{Images/FlashAdd3}}
\caption{Flash: comparison of unlimited solution with and without realizability limiter.}
\label{fig:flash3}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/Flash/}
\centering
\externaltikz{FlashAdd}{\input{Images/FlashAdd1}}
\caption{Flash: comparison of primitive and characteristic slope limiter.}
\label{fig:flash}
\end{figure}
%
\begin{figure}[h]
\def\localpath{./Images/Flash/}
\centering
\externaltikz{FlashAdd2}{\input{Images/FlashAdd2}}
\caption{Flash: comparison of primitive and characteristic slope limiter with realizability limiter.}
\label{fig:flash2}
\end{figure}



\clearpage
\subsection{Shadow}
\begin{itemize}
\item Domain: $[0,12]\times[0,6]$, $h=0.04$
\item Time: until steady state
\item Parameters: Let $O$ be the obstacle $O=[2,3]\times [0,2]$. Set $\sigma_s=0$; $\sigma_a=50$ on $O$, zero otherwise; source $q^{(0)}=0$; $q^{(1)} = 0$
\item Initial condition: $\psi^{(0)}=0$
\item Boundary conditions: top/bottom reflective, right vacuum, left Dirichlet with $\psi^{(0)}=1$, $\psi^{(1)}_x=0.99$, $\psi^{(1)}_y=0$
\end{itemize}

\begin{figure}[h]
\def\localpath{./Images/Shadow/}
\centering
\externaltikz{ShadowAdd3}{\input{Images/ShadowAdd3}}
\caption{Shadow: comparison of unlimited solution with and without realizability limiter.}
\label{fig:shadow3}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/Shadow/}
\centering
\externaltikz{ShadowAdd}{\input{Images/ShadowAdd1}}
\caption{Line source: comparison of primitive and characteristic slope limiter.}
\label{fig:shadow}
\end{figure}

\begin{figure}[h]
\def\localpath{./Images/Shadow/}
\centering
\externaltikz{ShadowAdd2}{\input{Images/ShadowAdd2}}
\caption{Line source: comparison of primitive and characteristic slope limiter with realizability limiter.}
\label{fig:shadow2}
\end{figure}

\clearpage

\subsection{Two Beams}
\begin{itemize}
\item Domain: $[0,7]\times[0,7]$, $h=0.05$
\item Time: 7
\item Parameters: $\sigma_a=\sigma_s=0$, $q^{(0)}=0$, $q^{(1)}=0$
\item Initial condition: $\psi^{(0)}=10^{-4}$, $\psi^{(1)}=0$
\item Boundary conditions: left $\psi^{(0)}=100$, $\psi^{(1)}_x=99.9$, $\psi^{(1)}_y=0$ in $0\times [3,4]$, bottom $\psi^{(0)}=100$, $\psi^{(1)}_x=0$, 
    $\psi^{(1)}_y=99.9$ in $[3,4]\times 0$,  else $\psi^{(0)}=10^{-4}$, $\psi^{(1)}=0$
\end{itemize}

\begin{figure}[h]
\def\localpath{./Images/TwoBeams/}
\centering
\externaltikz{TwoBeamsAdd3}{\input{Images/TwoBeamsAdd3}}
\caption{Two beams: comparison of unlimited solution with and without realizability limiter.}
\label{fig:Twobeams3}
\end{figure}
In Figure~\ref{fig:Twobeams3} we observe that in the absence of any limiting the numerical solution on the triangular unstructured mesh is unstable
and completely blows up. While the solution on the rectangular mesh does not completely blow up, it is highly oscillatory and loses realizability
on parts of the domain. Adding the realizability limiter stabilizes both numerical solutions; however, there are visible oscillations in 
the region where the beams collide.

\begin{figure}[h]
\def\localpath{./Images/TwoBeams/}
\centering
\externaltikz{TwoBeamsAdd}{\input{Images/TwoBeamsAdd1}}
\caption{Two beams: comparison of primitive and characteristic slope limiter.}
\label{fig:Twobeams}
\end{figure}
In Figure~\ref{fig:Twobeams} we note that limiting in the primitive or characteristic variables results in a solution that is not realizable on parts of the domain.
In the case of the unstructured triangular mesh the limiting in primitive variables results in an asymmetric solution. While
limiting in the characteristic variables leads to a better solution on both meshes, realizability cannot be guaranteed on the whole domain.
\begin{figure}[h]
\def\localpath{./Images/TwoBeams/}
\centering
\externaltikz{TwoBeamsAdd2}{\input{Images/TwoBeamsAdd2}}
\caption{Two beams: comparison of primitive and characteristic slope limiter with realizability limiter.}
\label{fig:Twobeams2}
\end{figure}

Finally in Figure~\ref{fig:Twobeams2}, the combination of limiting in characteristic variables plus
the realizability limiter leads to a superior solution on both meshes. However, we point out that in this test case
on the unstructured mesh, the numerical solution does not maintain a perfect symmetry when the beams collide. We have observed
that this is in part due to the fact that the transformations to and from the characteristic variables are ill conditioned
at the point when the beams collide.

dg.tex
\section{Runge-Kutta Discontinuous Galerkin Method}
\label{sec:DG}
%

\subsection{Spatial Discretization}
In this section, we describe the Runge-Kutta Discontinuous Galerkin method (RKDG) for solving~\eqref{eq:standard_form}--\eqref{eq:boundary_condition}. 
Following the approach outlined in a series of papers by Cockburn and Shu
\cite{CockburnShuIV,CockburnShuPk,CockburnShuP1,CockburnShuV}, we discretize in space using piece-wise
polynomials of degree $k$, that are allowed to be discontinuous at the cell interface. The time discretization is performed by a strong stability preserving explicit Runge-Kutta scheme of order $k+1$ \cite{gottlieb_shu_tadmor}.

In the following, let $\mathcal{T}_h$ be a partition of a polygonal computational domain $X$ and $K$ be
an element in $\mathcal{T}_h$ with boundary edges $e$. For each $t \in [0,T]$, we seek an approximate solution 
$U_h(x,t)$ in the DG space 
\[
V_h^k = \{v \in L^\infty(X): v|_{K} \in P^k(K), \forall K \in \mathcal{T}_h \}
\]
where $P^{k}(K)$ is the set of polynomials of degree $k$. We follow the Galerkin approach: first we multiply (\ref{eq:standard_form}) 
by a smooth test function $v$ and integrate by parts over each element. We replace the exact solution $U$ and smooth test function $v$ 
by the approximation solution $U_h$ and $v_h$ (both in $V_h^k$) respectively to obtain:
\begin{gather}
\label{eq:dweakform1}
\dfrac{d}{dt} \int_{K} U_h(x,t)v_h(x)\,dx + \sum_{e \in \partial K} \int_{e}\sF(U_h(x,t))\cdot
\bfn_{e,K} v_h(x)\,d\Gamma \nonumber\\ -\int_{K} \sF(U_h(x,t))	 \cdot \nabla v_h(x)\,dx 
= \int_{K} S\big(U_h(x,t) \big)v_h(x)\,dx \quad\forall v_h \in V_h^k,
 \label{eq:dweakform1a}
 \end{gather}
 where $\bfn_{e,K}$ is the outward unit normal to the boundary of $K$. We take $U_h(x,0)$ on each element to be the
 $L_2$-projection of the initial condition on $V_h^k$, i.e., 
 \begin{gather}
	 \int_{K} U_h(x,0)v_h(x)\,dx = \int_{K}u_0 v_h(x)\,dx, \quad \forall v_h \in V^k_h.
	 \label{eq:initial_condition}
 \end{gather}
The choice of a discontinuous basis implies that our approximate solution $U_h(x,t)$ is discontinuous across edges. 
In this case the normal trace $\displaystyle \sF(U_h(x,t))\cdot \bfn_{e,K}$ is not defined on the element boundary $\partial K$. 
We replace the normal trace by a numerical flux function $H_{e,K}(x,t)$ that depends
on the approximate solution from the interior and exterior of the element $K$. Thus we define
\begin{align}
\displaystyle H_{e,K}(x,t) = H_{e,K}(U_h(x^{int\{K\}},t), U_h(x^{ext\{K\}},t)), 
\label{eq:numerical_flux}
\end{align}
where
\begin{align*}
	U_h(x^{int\{K\}},t) &= \lim_{\xi\rightarrow (x,y) \in K}U_h(\xi,t),
\end{align*}
for the approximate solution defined from the interior of the element $K$ and
\begin{align*}
U_h(x^{ext\{K\}},t) = \left \{ \begin{array}{l} \gamma_h(x,t) \quad \textrm{ if } x \in \partial X,\\
                    \displaystyle\lim_{\xi\rightarrow (x,y) \notin K } U_h(\xi,t), \quad\textrm{otherwise} 
                    \end{array}
                    \right.
\end{align*}
for the points on the exterior of $K$.
Using the numerical flux \eqref{eq:numerical_flux}, the discrete weak formulation (\ref{eq:dweakform1})--(\ref{eq:initial_condition}) becomes
\begin{align}
\label{eq:dweakform2}
\dfrac{d}{dt} \int_{K} U_h(t,x)v_h(x)\,dx &+ \sum_{e \in \partial K}\int_{e} H_{e,K}(U_h(x^{int\{K\}},t),U_h(x^{ext\{K\}},t)) v_h(x)\,d\Gamma \nonumber\\
-\int_{K} \sF(U_h(x,t))\cdot \nabla v_h(x)\,dx                                                   
&= \int_{K} S(U_h(t,x))v_h(x)\,dx, \quad \forall v_h \in V^k_h,\\
\label{eq:dweakform2a}
\int_{K} U_h(x,0)v_h(x)\,dx &= \int_{K}u_0 v_h(x)\,dx, \quad \forall v_h \in V^k_h.
\end{align}
We choose the global Lax-Friedrichs flux
\begin{align}
\label{eq:LFFlux}
	H_{e,K}(a,b) = \dfrac{1}{2} \bigg[\sF(a)\cdot \bfn_{e,K} + \sF(b)\cdot \bfn_{e,K} - \alpha (b-a)\bigg],
\end{align}
where the numerical viscosity constant $\alpha$ is taken as the global estimate of the absolute value of the largest eigenvalue of the Jacobian \eqref{eq:jacobian}. 
For the $M_1$ model, we can take $\alpha=1$.

Boundary conditions \eqref{eq:boundary_condition} have to be incorporated via the quantities $U_h(x^{ext\{K\}},t)$ in the edge integral, if the edge $e$ is part of the boundary. We take the simplest approach possible by adding so-called ``ghost cells'' which are then filled with the corresponding values for $U$. 
Note, however, that the validity of this approach, due to its inconsistency with 
the original boundary conditions  \eqref{eq:boundary_condition}, is not entirely non-controversial, but the question of appropriate boundary conditions for
moment models is an open problem \cite{pomraning1964variational,Larsen1991,Rulko1991,Struchtrup2000,levermore2009boundary} which is not explored here. 

Equations (\ref{eq:dweakform2})--(\ref{eq:dweakform2a}) can be written as a system of ODEs on each element after inverting the mass matrix in~\eqref{eq:dweakform2}.
Indeed, let $\displaystyle \{\varphi_1, \varphi_2, \cdots ,\varphi_{N_k} \}$ denote a basis of the space of 
polynomials of degree $k$ on cell $K$, where $\displaystyle N_k = \dfrac{1}{2}(k+1)(k+2)$. On each element, the
DG approximate solution $U_h$ of the components has the form
\begin{align}
\label{eqn:solution_form}
\psi^{(0)} = \sum_{i=1}^{N_k}\alpha^{\psi^{(0)}}_i \varphi_i,\quad \psi_x^{(1)} = \sum_{i=1}^{N_k}\alpha^{\psi^{(1)}_x}_i \varphi_i, \textrm{ and }
\psi_y^{(1)} = \sum_{i=1}^{N_k}\alpha^{\psi^{(1)}_y}_i \varphi_i,
\end{align}
where $\alpha^{\psi^{(0)}}_i, \alpha^{\psi^{(1)}}_x, \alpha^{\psi^{(1)}}_y$ are unknowns to be determined.
Using the form of the approximate solution in (\ref{eqn:solution_form}), we can write (\ref{eq:dweakform2})-(\ref{eq:dweakform2a}) in matrix form:
\begin{gather}
\label{eqn:weakform_m}
\bfM\dfrac{\partial}{\partial t} \alpha_h + \bfH \alpha_h - \bfF\alpha_h = \bfS \alpha_h \\
\label{eqn:initial_condition}
\bfM \alpha^0_h = \bfU_0
\end{gather}
where $\alpha_h$ is the vector of solution coefficients and 
\begin{gather}
	(\bfM)_{ij} = \int_{K}\varphi_j \varphi_i,\quad (\bfH)_{i} = \sum_{e \in \partial K} \int_{e}H_{e,K}(x,t)\varphi_i,
    \quad (\bfF)_{i} = \int_{K}\mathcal{F}(U_h)\cdot \nabla \varphi_i\\
	(\bfS)_i = \int_{K} S(U_h)\varphi_i, \textrm{ and } (\bfU_0)_i = \int_{K}u_0\varphi_i.
\end{gather}
The complete coefficient vector is given by $\displaystyle \alpha_h = \{\alpha^{\psiz}_1, \cdots, \alpha^{\psiz}_{N_k}, \alpha^{\psio_x}_{1},\cdots,\alpha^{\psio_x}_{N_k},\alpha^{\psio_y}_{1},\cdots,\alpha^{\psio_y}_{N_k} \}$, where $\displaystyle \alpha^{\psiz}_i, \alpha^{\psio_x}_i, \alpha^{\psio_y}_i$ 
are the coefficients of the numerical approximation to the zeroth and first order moments, respectively. 
We can write \eqref{eqn:weakform_m}--\eqref{eqn:initial_condition} in the form
\begin{gather}
\label{eq:ode1}
\dfrac{d}{dt} \alpha_h = \mathcal{L}_h(\alpha_h), \quad \textrm{on } X \times (0,T) \\
\label{eq:ode2}
\alpha_h(x,0) = \alpha^{0}_h 
\end{gather}
where
\begin{gather}
\label{eq:ode3}
\mathcal{L}_h(\alpha_h) = \bfM^{-1}(\bfS\alpha_h + \bfF \alpha_h - \bfH \alpha_h).
\end{gather}
We approximate the solution by discontinuous quadratic polynomials in space and the third-order strong stability preserving Runge-Kutta time discretization scheme 
proposed in \cite{CWShu,Shu1988}, also known as the Shu-Osher scheme. Let $\{t^n\}_{n=0}^{N}$ be a partition of $[0,T]$ and let $\Delta t = t^{n+1} - t^n$, $n = 0, \cdots, N-1$, then the time
stepping scheme updating the coefficients of the DG polynomials can be written as:
\begin{center}
	\begin{itemize}
        \item  Set  $\alpha^0 =\bfM^{-1}\bfU_0$;
			\item  For  $n = 1, \cdots, N-1$  compute  $\alpha_h^{n+1}$ as follows:
				\begin{enumerate}
					\item $\displaystyle \alpha_h^{(1)} = \alpha_h^n + \Delta t_{n} \mathcal{L}_h(\alpha^n_h)$
					\item $\displaystyle \alpha_h^{(2)} = \dfrac{3}{4}\alpha_h^n + \dfrac{1}{4}(\alpha_h^{(1)} + \Delta t_{n} \mathcal{L}_h(\alpha_h^{(1)}))$
					\item $\displaystyle \alpha_h^{(3)}= \dfrac{1}{3}\alpha_h^{n} + \dfrac{2}{3}(\alpha_h^{(2)} + \Delta t_{n} \mathcal{L}_h(\alpha_h^{(2)}))$
					\item Set $\alpha_h^{n+1} = \alpha_h^{(3)}$,
				\end{enumerate}
		\end{itemize}
\end{center}
This Runge-Kutta method is a convex combination of (iterated) forward Euler steps. Using the convexity of the realizability domain, one can show 
that it preserves realizability under a specific CFL-condition. We use this property to achieve high order also in time without the need of dealing with complicated time-discretizations in the proof of realizability preservation.

\subsection{Quadrature Rules}
%
The assembly of the discrete operator $\displaystyle \mathcal{L}_h(\alpha_h)$ (\ref{eq:ode1}) is done using numerical
quadrature that is exact for polynomials of degree $2k+1$ for edge integrals, and degree $2k$ for volume
integrals, respectively for both the rectangular and triangular meshes. In addition, we will need a modified Gaussian quadrature 
rule to construct the realizability limiter. This rule consists of quadrature points in the interior of the cell and on its boundary where the latter form one-dimensional quadratures on the element 
edges. This is crucial to balance the different appearing types of spatial integrals in the proof of the realizability-preserving property.
These rules are described in the following section that closely follows the construction in \cite{ZhangShu12}.

\subsubsection{Triangles}
%
We start with two quadrature rules on the interval $[-\frac{1}{2},\frac{1}{2}]$, a Gaussian rule
$$
v^\beta = \frac{1}{2}\begin{pmatrix}
-\sqrt{\frac{3}{5}}&0&\sqrt{\frac{3}{5}}
\end{pmatrix},\quad w^\beta = \frac{1}{18} \begin{pmatrix}
5 & 8 & 5
\end{pmatrix}
$$
and a Gauss-Lobatto rule
\begin{align}
\label{eq:GLrule}
u^\alpha = \frac{1}{2}\begin{pmatrix}
-1&0&1
\end{pmatrix},\quad
\hat{w}^\alpha = \frac{1}{6} \begin{pmatrix}
1 & 4 & 1
\end{pmatrix}.
\end{align}
We map the tensor product rule of those Gauss and Gauss-Lobatto quadrature rules from the square $[-\frac{1}{2},\frac{1}{2}]^2$ onto the triangle $K$ with vertices $V_1^K,\ldots,V_3^K$ using the three degenerate projections 
\begin{align*}
g_1(u,v) &= (\frac{1}{2}+v)V^K_1 + (\frac{1}{2}+u)(\frac{1}{2}-v)V^K_2+(\frac{1}{2}-u)(\frac{1}{2}-v)V^K_3\\
g_2(u,v) &= (\frac{1}{2}+v)V^K_2 + (\frac{1}{2}+u)(\frac{1}{2}-v)V^K_3+(\frac{1}{2}-u)(\frac{1}{2}-v)V^K_1\\
g_3(u,v) &= (\frac{1}{2}+v)V^K_3 + (\frac{1}{2}+u)(\frac{1}{2}-v)V^K_1+(\frac{1}{2}-u)(\frac{1}{2}-v)V^K_2.
\end{align*}
Let $\bar{U}_K$ denote the cell average of the numerical solution $U_h$ represented by the DG polynomial $p_K$ on triangle $K$, following \cite{ZhangShu12} it holds that 
\begin{align*}
\bar{U}_K = \frac{2}{3}\sum\limits_{i=1}^{3}\sum\limits_{\alpha=1}^{3}\sum\limits_{\beta=1}^{3}p_K\left(g_i\left(\hat{u}^\alpha,v^\beta\right)\right)\left(\frac{1}{2}-v^\beta\right)w^\alpha\hat{w}^\beta.
\end{align*}
Note that the Gauss-Lobatto points on the edge are always taken twice (see Figure \ref{fig:Quadrature}).
%
\begin{figure}[htbp]
\centering
\settikzlabel{fig:quadrulestrianglePos}
\settikzlabel{fig:quadrulestriangleWeights}
\externaltikz{quadrulestriangle}{\input{Images/quadrulestriangle}}  
\caption{Position of the quadrature nodes under the three projections $g_i$. The color corresponds to the weight $\tilde{w}^\gamma$ at the quadrature node.}
\label{fig:Quadrature}
\end{figure}

The authors of \cite{ZhangShu12} showed that the cell mean can be constructed as the following convex combination of inner and boundary points:
\begin{align}
\label{eq:Cellmean}
\bar{U}_{K} = \sum\limits_{i=1}^{3}\frac{2}{3}w^\beta\hat{w}^{\alpha=1} U_{i,\beta}^{int}+\sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma U_\gamma^{inner}
\end{align}
where $U_{\gamma}^{inner}$ and $U_{i,\beta}^{int}$ represent the evaluation of the DG polynomial $p_K$ at the respective interior and boundary quadrature nodes 
$(x_{\gamma}^{inner},y_{\gamma}^{inner})$ and $(x_{i,\beta}^{int},y_{i,\beta}^{int})$ (compare Figure~\ref{fig:Quadrature}), and $L = 3(N-2)(k+1) = 9$ is the number of inner points. With this we obtain a quadrature rule with $18$ points which is accurate for polynomials of order $2k-1=3$.

The flux term in \eqref{eq:dweakform2} should be approximated by the $1D$ $(k+1)$ point Gauss quadrature with weights $w^\beta$:
\begin{align}
\label{eq:FluxGauss}
\sum\limits_{e\in \partial K}\int_e H_{e,K}(x,t) \varphi_i^K(x)~d\Gamma = \sum\limits_{i=1}^3\sum\limits_{\beta = 1}^{3} H_{e_i,K}(U_{i,\beta}^{int},U_{i,\beta}^{ext},t) \varphi_i^K(x_{i,\beta}^{int},y_{i,\beta}^{int})w^\beta l^i_K
\end{align}
where $l^i_K$ is the length of edge $e_i^K$. 

\subsubsection{Rectangles}
%
For the assembly of the discrete operator in the rectangular DG code we use a tensor quadrature rule of the 1D  $4$-point Gauss-Lobatto rule on $[-\frac12,\frac12]$:
$$
u^\alpha = \cfrac{1}{2\sqrt{5}}\begin{pmatrix}
-\sqrt{5}&-1&1&\sqrt{5}
\end{pmatrix},\quad
\hat{w}^\alpha = \frac{1}{12} \begin{pmatrix}
1 & 5 & 5 & 1
\end{pmatrix}.
$$
However for the proof of our main theorem we only need the three point rule \eqref{eq:GLrule} which will give a weaker CFL-condition.

\subsection{Slope Limiting}
In order to enforce stability \cite{Cockburn,hesthaven2007nodal} and to mitigate numerical oscillations caused by the Gibbs phenomenon we apply a slope limiter
to every stage of the Runge-Kutta time stepping scheme. We implement the slope limiter detailed 
by Cockburn and Shu in \cite{CockburnShuIV}
\begin{gather}
	\overline{m}(a_1,a_2,\cdots, a_n) = \begin{cases} a_1 & \mbox{ if } |a_1| < M (\Delta x)^2 \\
										m(a_1,a_2,\cdots,a_n) & \mbox{ otherwise,}
\end{cases}
\label{eq:dg_slope_limiter}
\end{gather}
where $M$ is an estimate of the second derivative of the solution, $\Delta x$ is the characteristic length of
each element and $\overline{m}$ is the TVB modified
minmod function. The limiter can be applied component-wise to the primitive variables
or in the characteristic variables. For limiting in the characteristic variables we first construct a matrix
$\mathcal{R}$ that diagonalizes the directional Jacobian in the normal direction $\bfn = (n_x,n_y)$ evaluated at the mean in each element
\begin{gather}
\label{eq:trafo}
	\mathcal{R}^{-1} (n_x J_x(\overline{U}_k) + n_y J_y(\overline{U}_k))\mathcal{R} =  \Lambda.
\end{gather}
Applying the limiter in the characteristic variables ensures that the solution is total variation bounded in the means (TVBM). 



introduction.tex
\section{Introduction}
%
The $M_1$ model of radiative transfer is a nonlinear system of hyperbolic balance laws, and reads
\begin{subequations}
\begin{align}
    \partial_t \psiz + \nabla_x\cdot \psio &= -\sigma_a \psiz + q^{(0)} \label{eq:m1_1} \\
    \partial_t \psio + \nabla_x \cdot \psit(\psiz,\psio) &= -(\sigma_s+\sigma_a) \psio + q^{(1)}.\label{eq:m1_2}
\end{align}
\end{subequations}

The quantities $\psiz$, $\psio$, $\psit$ are, respectively, the zeroth (particle density), first (mean velocity) and second moment (pressure) over the unit sphere $S^2$ of the angular flux $\psi$, 
$$
\psiz := \int_{S^2} \psi(\Omega)d\Omega,\quad \psio := \int_{S^2} \Omega\psi(\Omega)d\Omega,\quad\psit := \int_{S^2} \Omega\Omega^T\psi(\Omega)d\Omega, 
$$
and are, respectively, a scalar, a vector, and a matrix. The system is supplemented by the closure condition
\begin{equation}
\label{eq:closure}
\psit(\psiz,\psio) = D(\frac{\psio}{\psiz}) \psiz,
\end{equation}
where 
\begin{gather}
D(n) = \frac{1-\chi(|n|)}{2}\text{id} + \frac{3\chi(|n|)-1}{2}\frac{nn^T}{|n|^2}\text{ and }\chi(f) = \frac{3+4f^2}{5+2\sqrt{4-3f^2}}~\text{ for } f\in[0,1].
\label{eq:D_exp}
\end{gather}
Here, $D$ is a matrix, $n=\frac{\psio}{\psiz}$ is a vector, and $f=|n|$ is a scalar. The quantity $\chi(f)$ is called the Eddington factor \cite{Lev84,BruHol01}. Let $X \subset \mathbb{R}^2$ be a bounded polygonal domain, then the system \eqref{eq:m1_1} -\eqref{eq:m1_2} 
can be written as a general first-order system of balance laws
\begin{gather}
	\dfrac{\partial U}{\partial t} + \nabla \cdot \sF   =  S,\text{ in } X \times (0,T) \label{eq:standard_form}\\
U(x,y,0) = u_{0}(x,y),\text{ for } (x,y)  \in X \label{eq:initial_condition}\\
\end{gather}
where $\displaystyle \sF = [F,G]$ and
$$
U = \begin{bmatrix} \psiz \\ \psio_x \\ \psio_y \end{bmatrix},\ 
F = \begin{bmatrix} \psio_x \\ \psit_{xx} \\ \psit_{xy} \end{bmatrix},\ 
G = \begin{bmatrix} \psio_y \\ \psit_{yx} \\ \psit_{yy} \end{bmatrix},\ 
S = \begin{bmatrix} -\sigma_a \psiz + q^{(0)} \\ -(\sigma_s+\sigma_a) \psio_x + q^{(1)}_x \\ -(\sigma_s+\sigma_a) \psio_y + q^{(1)}_y \end{bmatrix},
$$
keeping in mind the closure relation \eqref{eq:closure}. Here, $\psio_x$ and $\psio_y$ denote the first and second component of $\psio$, respectively, similarly for $\psit_{xx}$ to $\psit_{yy}$. 

The system has to be complemented with boundary conditions of the form
\begin{align}
U(x,y,t) = \gamma(x,y,t),\text{ for } (x,y) \in I(\partial X,U) \times (0,T) \label{eq:boundary_condition},
\end{align}
where the operator $I$ returns the influx boundaries, i.e.\  those parts of the boundary where information is transported into the domain \cite{toro2009riemann}. Whether a part of the boundary is an influx boundary also depends on the solution itself, since the information direction can be read off the sign of the eigenvalues of the directional Jacobian in normal direction at these points. In our numerical experiments, we consider
compactly supported initial data which do not reach the boundary
before the final time. We prescribe Dirichlet boundary conditions
identical to the initial conditions, hence the numerical solution is
unaffected by the boundaries.

The expression for $D$ in \eqref{eq:D_exp} comes from closing the moment system by an entropy closure using the entropy for photons. See for example \cite{Dre87,DubFeu99,Frank07,jaynes1957info,Lev84,levermore1996moment,MullerRuggeri,Ore55,Ros54} and references therein for more information. These references also discuss many properties which make the $M_1$ model, and entropy closures in general, quite appealing. Among these are symmetric hyperbolicity, i.e.\  the system can be transformed into a symmetric hyperbolic system, and a natural entropy-entropy flux pair, both of which ensure some level of well-posedness \cite{toro2009riemann,LeVeque2002}. The hyperbolicity is of main interest in this work, and will be discussed in detail below. We especially focus on its connection to realizability -- the fact that a moment vector can be reproduced by a non-negative particle density.

When the $M_1$ model is discretized using a monotone first-order scheme for hyperbolic equations (e.g. the Lax-Friedrichs scheme), it can be shown that the numerical solution will never leave the set of realizable moments \cite{Olbrant12}. This means that starting with a realizable initial condition, the numerical solution will be realizable at every time step. This is in particular required, since the model is not well-defined outside of the realizability set.

Unfortunately, and this is the main topic of this paper, higher-order numerical schemes do not automatically preserve this property, as will be discussed later. We use the Discontinuous Galerkin (DG) method, which provides a general framework to construct numerical schemes of arbitrary approximation order to solve hyperbolic balance laws. The original DG method was introduced in $1973$ by Reed and Hill \cite{reedhill} for neutron transport and has been developed further to the Runge-Kutta Discontinuous Galerkin (RKDG) method by Cockburn et al.\ in a series of papers \cite{CockburnShuIV,CockburnShuPk,CockburnShuP1,CockburnLinShu,CockburnShuV}. A $(k+1)st$ order RKDG method uses a piecewise-polynomial approximation of degree $k$ in space and a ($k+1)st$ order strong stability preserving explicit Runge-Kutta scheme in time. 

A DG scheme has to be supplemented by a limiting strategy. This typically consists of a slope limiter which ensures stability of the solution. A standard example of such a limiter is the TVBM limiter \cite{CockburnLinShu}. For a system of balance laws, limiting the
conserved variables component-wise can result in oscillations, due to the Gibbs phenomenon \cite{Shu1998}. Instead, the limiter has to be applied in characteristic variables, i.e.\ \ Riemann invariants that are obtained by diagonalizing the system Jacobian \cite{biswas1994parallel,CockburnShuV}. In general, limiting in the characteristic variables gives superior results \cite{review_article,Olbrant12}. However, there are numerous examples
for which the numerical solution obtained by limiting in the characteristic variables and component-wise limiting are comparable, e.g. \cite{biswas1994parallel,li2011central,liu2007runge,liangdiscontinuous,liu_xu_gas_kinetic}. In these situations, since the transformation to and from the characteristics requires additional computational effort, component-wise 
limiting is faster. 

It has been shown in \cite{Olbrant12} that a robust numerical approximation of the $M_1$
model in 1D requires an additional limiter: a realizability limiter
which ensures that the model remains well-posed. This has been extended to the general, one-dimensional, $M_N$ case in \cite{Schneider2015a}. In the current work, we extend this realizability limiting strategy to two dimensions. In addition, we demonstrate that in the absence of the realizability limiter, the TVBM limiter applied in the characteristic variables
leads to a qualitatively superior solution compared to component-wise limiting; however, this limiter alone is not enough to guarantee moment realizability of the scheme.

The $M_1$ model therefore might serve as a good benchmark test for DG implementations, because it absolutely requires very careful limiting. 
We arrive at this conclusion experimentally, by running several well-known radiation transport test cases with two independently developed third-order 
DG schemes on unstructured triangular grids and regular rectangle grids. 

The rest of this paper is organized as follows. 
In Section~\ref{sec:M1}, we describe in detail the concept of realizability, and how it is connected to the well-posedness of the $M_1$ model. 
The DG implementations are described in Section~\ref{sec:DG}. The realizability limiter will be developed and analyzed in 
Section \ref{sec:limiter}. Section \ref{sec:numerical_results} contains the results from the test cases.

numerics.tex
\section{Numerical Results}
\label{sec:numerical_results}
In this section, we compare the performance of the RKDG schemes on unstructured triangular meshes
and uniform rectangular meshes for various combinations of limiters. Specifically we consider the limiter combinations given in Table~\ref{tab:LimiterCombinations}.
The time step is taken as indicated by~\eqref{eqn:cfl}.
\begin{table}[h]
\centering
\begin{tabular}{l l}
Abbr. & Explanation\\
\midrule
SL$M$$\square/\triangle$ & Slope limiter in primitive variables with the constant $M$ specified\\
CL$M$$\square/\triangle$ & Slope limiter in characteristic variables with the constant $M$ specified\\
SRL$M$$\square/\triangle$ & Slope limiter in primitive variables + realizability limiter\\
CRL$M$$\square/\triangle$ & Slope limiter in characteristic variables + realizability limiter\\
\end{tabular}
\caption{Abbreviation of limiter combinations. Setting $M=\infty$ is equivalent to disabling the slope limiter. $\square$ and $\triangle$ correspond to the rectangular and triangular meshes, respectively.}
\label{tab:LimiterCombinations}
\end{table}
To ensure that our research is reproducible, the codes used in this publication can be found on \textit{GitHub}, see \cite{Chidyagwai2016}.

In the numerical results that follow we denote the mesh size on both unstructured triangular meshes
and uniform rectangular meshes by $h$. All plots show either the zeroth moment $\psi^{(0)}$ or the norm of the first normalized moment $\vert\phi^{(1)}\vert = \vert\frac{\psi^{(1)}}{\psi^{(0)}}\vert$. The loss of realizability is indicated in the latter by the use of black ($\vert\phi^{(1)}\vert>1$) and white ($\psi^{(0)}<0$) colors. 

Before we describe the test cases and the results, we detail some of the steps in the algorithm, and the computation of the errors.


{\bfseries Fixing a DG solution which is not realizable.}
In test cases where the realizability limiter is not used (SL$M$/CL$M$) the solution may violate the realizability condition and thus the flux function cannot 
be evaluated. We therefore compute a modified numerical solution which we use only to compute the flux, i.e.\ the actual solution in each cell remains the same. 
We modify the values of the numerical solution locally at each Gauss point
as follows: let $\epsilon$ be a specified numerical tolerance, if $\psi^{(0)} < \epsilon$ then 
we set $\psi^{(0)} = \epsilon$. In addition, if the normalized first moment 
$$
f :=  |\frac{\psi^{(1)}}{\psi^{(0)}}| > (1 - \epsilon),
$$ 
we modify the values of first moment as follows: $\displaystyle \psi^{(1)}_x = \dfrac{\psi^{1}_x}{f}(1-\epsilon)$ and 
$\displaystyle \psi^{(1)}_y = \dfrac{\psi^{1}_y}{f}(1-\epsilon)$. We choose a tolerance of $\epsilon = 10^{-12}$.

{\bfseries Transformation to characteristic variables.}

In the case of the triangular meshes the matrices $\mathcal{R}$ and $\mathcal{R}^{-1}$ from \eqref{eq:trafo} are computed using the 
GSL linear algebra package \cite{Gough_GSL} while they can be precomputed analytically in case of the rectangular grid. 
The transformation to and from the characteristic variables is achieved by left multiplying the variables by $\mathcal{R}^{-1}$ and $\mathcal{R}$, respectively. 

In case that we need the transformation to characteristic variables and do not use the realizability limiter (i.e. CL$M$), we use the same realizability fix as in the flux function, but for the computation of the transformation matrices only. This ensures that transformation matrices are invertible due to the strict hyperbolicity of the moment system in the interior of the realizable set. The transformation matrix has a condition number of roughly $2/\epsilon$, where $\epsilon$ is the distance of the (cell mean of the) normalized first moment to the realizability boundary. Note that this implies that close to the realizability boundary, significant round-off errors may occur.
However, we have performed numerical experiments with the linesource test case with $\epsilon=10^{-5}, 10^{-12}, 10^{-14}$, and observed no significant effect of the cutoff parameter.

{\bfseries Comparison to a reference solution.}
We note that for a general $f>0$ we have that $$\lVert\log_{10}(f)\rVert_{H^1(K)} = \lVert\log_{10}(f)\rVert_{L^2(K)}+\lVert\frac{1}{\log_{10}f}\nabla_xf\rVert_{L^2(K)}.$$ 
We can therefore use the evaluation of the DG polynomial and its first-order derivatives in every element $K$ to calculate its logarithmic Sobolev representation locally.
 
The reference solution is computed using a first-order Lax-Friedrichs finite-volume scheme on an equidistant rectangular grid with $8192\times 8192$ grid points. The gradients for the reference solution are obtained using centered finite difference formulas. Finally, the argument in the norm is evaluated on every point of the reference grid and integrated using a quadrature with equal weights at exactly those reference points (rectangle rule). 


{\bfseries Choosing the limiter parameter $M$.}
The solutions of the benchmark test cases that follow are plotted on a logarithmic scale. We thus choose the value of $M$ in \eqref{eq:dg_slope_limiter}, using the logarithmic 
Sobolev norm error, as
\begin{align}
\label{eq:logSob}
M = \operatorname{arg min}\limits_{M\geq 0} \lVert \log_{10}\left(\psi^{(0)}_{M}\right)-\log_{10}\left(\psi^{(0)}_{ref}\right)\rVert_{H^1(X)}.
\end{align}
Since every evaluation of \eqref{eq:logSob} requires to solve the full system of equations, finding the true minimum using e.g. a gradient method is not feasible. We therefore restricted ourself to the discrete set of (almost) logarithmically spaced values $M\in\mathcal{M}:=\{0.1,0.2,0.5,1,2,10,22,46,100,150\}$.

\subsection{Failure of the realizability limiter}
%
We want to extend the investigation in \cite{Schneider2015a} to our two-dimensional setup. 
Define the two moment vectors 
\begin{align*}
U_0 &= (1-\xi)\begin{pmatrix}
1\\1\\0
\end{pmatrix} + \xi\begin{pmatrix}
1\\0\\0
\end{pmatrix},&
U_1 = 10^{-6}\cdot\left((1-\xi)\begin{pmatrix}
1\\0\\1
\end{pmatrix} + \xi\begin{pmatrix}
1\\0\\0
\end{pmatrix}\right).
\end{align*}
Both of them satisfy $\lvert \frac{\psio}{\psiz}\rvert = 1-\xi$, i.e. the distance to the realizability boundary can be controlled by the parameter $\xi$. Due to the convexity of the realizable set every convex combination of $U_0$ and $U_1$ will be realizable as well. We can thus define the vector of functions
\begin{align}
\label{eq:u-limit-test}
U(x,y) &= (1-\lambda(x,y))U_0 + \lambda(x,y)U_1,\\
\lambda(x,y) &= \cfrac{\cos\left(2\left(x+y\right)\pi\right)+1}{2}\in[0,1],\nonumber
\end{align}
which is realizable for all $(x,y)\in[0,1]^2$.

Now we discretize the domain, for simplicity, using the structured rectangular mesh with equidistant step sizes. We project $U(x,y)$ onto the DG basis of degree $k$, $k\in\{0,1,2\}$ and apply the realizability-preserving limiter from Section~\ref{sec:limiter}.
Then we refine the grid and calculate $L_1$- and $L_\infty$-errors for the zeroth moment $\psiz$ as 
\begin{align*}
E^1_h = \int_{0}^{1}\int_{0}^{1} \abs{\psiz_a(x,y)-\psiz_h(x,y)}~dxdy, \quad E^\infty_h = \max\limits_{(x,y)\in[0,1]^2}\abs{\psiz_a(x,y)-\psiz_h(x,y)},
\end{align*}
where $\psiz_a$ denotes the exact zeroth moment of $U$ and $\psiz_h$ its limited DG polynomial.
In practice, the integral in $E^1_h$ is approximated by a high-order tensor Gaussian quadrature on every cell and the $\max$ in $E^\infty_h$ is taken over the same quadrature nodes.

The observed convergence order $\nu$ is defined by
\begin{equation}
 \frac{E^p_{h1}}{E^p_{h2}} = \left( \frac{h_1}{h_2} \right)^\nu
\label{eq:conv-order}
\end{equation}
where for $i \in \{1, 2\}$, $E^p_{hi}$ is the error $E^p_h$ for the
numerical solution using cell size $h_i$, for $p \in \{1, \infty\}$.

We found that taking $\xi \in [0, 10^{-4}]$ places the moment curve
$U(x,y)$ close enough to the boundary of realizability that the realizability
limiter was active for every number of cells we considered.
In Table~\ref{tab:limiter-test-4} we show convergence rates for 
$\xi=10^{-4}$. These results show the expected convergence order $k+1$ in all cases.
In this table we include the column $\theta_{\max}$, which gives the
maximum value of $\theta$ from the realizability limiter over all spatial
cells.
That $\theta_{\max}$ is nonzero in each row indicates that the
realizability limiter was active for every reconstruction.

\begin{table}
\centering
\begin{tabular}{r r@{.}l c r@{.}l c r@{.}l r@{.}l c r@{.}l c r@{.}l}

 & \multicolumn{8}{c}{$k=1$} & \multicolumn{8}{c}{$k=2$}\\
\cmidrule(r){2-9} \cmidrule(r){10-17}
$1/h$ & \multicolumn{2}{c}{$E^1_h$} & $\nu$
 & \multicolumn{2}{c}{$E^\infty_h$} & $\nu$
 & \multicolumn{2}{c}{$\theta_{\max}$}
 & \multicolumn{2}{c}{$E^1_h$} & $\nu$
 & \multicolumn{2}{c}{$E^\infty_h$} & $\nu$
 & \multicolumn{2}{c}{$\theta_{\max}$} \\ \midrule
 
 5 & 5 & 535e-02 & --- & 2 & 790e-01 & --- & 3 & 210e-01& 1 & 787e-02 & --- & 6 & 567e-02 & --- & 1 & 783e-01\\
10 & 1 & 252e-02 & 2.1 & 1 & 269e-01 & 1.1 & 3 & 960e-01& 1 & 483e-03 & 3.6 & 1 & 153e-02 & 2.5 & 5 & 305e-02\\
20 & 2 & 774e-03 & 2.2 & 3 & 831e-02 & 1.7 & 4 & 117e-01& 1 & 382e-04 & 3.4 & 1 & 495e-03 & 2.9 & 1 & 391e-02\\
40 & 6 & 500e-04 & 2.1 & 1 & 001e-02 & 1.9 & 4 & 154e-01& 1 & 551e-05 & 3.2 & 1 & 878e-04 & 3.0 & 3 & 519e-03\\
80 & 1 & 568e-04 & 2.1 & 2 & 531e-03 & 2.0 & 4 & 164e-01& 1 & 881e-06 & 3.0 & 2 & 350e-05 & 3.0 & 8 & 824e-04\\
160 & 3 & 847e-05 & 2.0 & 6 & 345e-04 & 2.0 & 4 & 166e-01& 2 & 332e-07 & 3.0 & 2 & 938e-06 & 3.0 & 2 & 205e-04\\
320 & 9 & 526e-06 & 2.0 & 1 & 587e-04 & 2.0 & 4 & 166e-01& 2 & 909e-08 & 3.0 & 3 & 673e-07 & 3.0 & 5 & 431e-05\\


\end{tabular}
\caption{$L^1$- and $L^\infty$-errors and observed convergence order $\nu$
for the first component of the realizability-limited, piece-wise linear and
quadratic reconstructions of $U(x,y)$ from \eqref{eq:u-limit-test} with
$\xi = 10^{-4}$.}
\label{tab:limiter-test-4}
\end{table}

Pushing $U$ closer to the realizability boundary ($\xi\leq 10^{-5}$) degrades the convergence order of the piecewise-linear reconstruction, verifying the mentioned convergence problems. This is demonstrated in Table~\ref{tab:limiter-test-0} for the extreme case $\xi = 0$, i.e. moments are placed on the realizability boundary. Surprisingly, the third-order approximation has the full convergence order for all tested $\xi\in[0,10^{-3}]$.
Note that this is different to the results in \cite{Schneider2015a} where the third-order approximation degrades to second order. The reason for this is still unclear but might be due to the high non-linearity of the $M_3$ model used in \cite{Schneider2015a}. We checked several variants of the test (by choosing different vectors on the unit sphere instead of $(1,0,0)^T$ and modifying the weight $10^{-6}$ in front of $U_1$) and all showed the same behavior. 

\begin{table}
\centering
\begin{tabular}{r r@{.}l c r@{.}l c r@{.}l r@{.}l c r@{.}l c r@{.}l}

 & \multicolumn{8}{c}{$k=1$} & \multicolumn{8}{c}{$k=2$}\\
\cmidrule(r){2-9} \cmidrule(r){10-17}
$1/h$ & \multicolumn{2}{c}{$E^1_h$} & $\nu$
 & \multicolumn{2}{c}{$E^\infty_h$} & $\nu$
 & \multicolumn{2}{c}{$\theta_{\max}$}
 & \multicolumn{2}{c}{$E^1_h$} & $\nu$
 & \multicolumn{2}{c}{$E^\infty_h$} & $\nu$
 & \multicolumn{2}{c}{$\theta_{\max}$} \\ \midrule
 
 5 & 6 & 716e-02 & --- & 3 & 594e-01 & --- & 5 & 603e-01& 1 & 787e-02 & --- & 6 & 567e-02 & --- & 1 & 783e-01\\
10 & 1 & 532e-02 & 2.1 & 1 & 269e-01 & 1.5 & 3 & 960e-01& 2 & 033e-03 & 3.1 & 1 & 153e-02 & 2.5 & 5 & 305e-02\\
20 & 3 & 143e-03 & 2.3 & 3 & 831e-02 & 1.7 & 4 & 117e-01& 1 & 577e-04 & 3.7 & 1 & 495e-03 & 2.9 & 1 & 391e-02\\
40 & 6 & 966e-04 & 2.2 & 1 & 001e-02 & 1.9 & 4 & 154e-01& 1 & 614e-05 & 3.3 & 1 & 878e-04 & 3.0 & 3 & 519e-03\\
80 & 1 & 626e-04 & 2.1 & 2 & 531e-03 & 2.0 & 4 & 164e-01& 1 & 901e-06 & 3.1 & 2 & 350e-05 & 3.0 & 8 & 829e-04\\
160 & 3 & 920e-05 & 2.1 & 6 & 345e-04 & 2.0 & 4 & 166e-01& 2 & 338e-07 & 3.0 & 2 & 938e-06 & 3.0 & 2 & 217e-04\\
320 & 9 & 625e-06 & 2.0 & 3 & 008e-04 & 1.1 & 4 & 167e-01& 2 & 912e-08 & 3.0 & 3 & 673e-07 & 3.0 & 8 & 294e-05\\
\end{tabular}
\caption{$L^1$- and $L^\infty$-errors and observed convergence order $\nu$
for the first component of the realizability-limited, piece-wise linear and
quadratic reconstructions of $U(x,y)$ from \eqref{eq:u-limit-test} with
$\xi = 0$.}
\label{tab:limiter-test-0}
\end{table}

Full convergence tests for this scheme in one and two dimensions can be found in \cite{Schneider2015a,Schneider2016}.

%%%%%%%%%%%%%%%%%%%%%%%%
%Linesource
\label{sec:results}
\subsection{Line Source}
%\comment{M = 22}
Our first benchmark is the line-source test as proposed in \cite{Brunner2005b}. It is a Green's function problem, where a pulse of particles is emitted from a line in an infinite medium. This corresponds to an initial condition of the form $\psi^{(0)}(x,y) = c\cdot\delta(x,y)$, whose lack of regularity makes this problem especially hard to tackle for the method of moments, including the $M_1$ model, see e.g. \cite{Garrett2014}. The parameters of this test are as follows:
\begin{itemize}
\item Domain: $[-0.5,0.5]^2$ with $h =0.004$
\item Time $T=0.45$
\item Parameters: $\sigma_a=\sigma_s=0$, $q^{(0)}=0$, $q^{(1)}=0$
\item Initial condition: smoothed Dirac $\psi^{(0)}(x,y) = \max(\exp(-10\frac{x^2+y^2}{\sigma^2}),10^{-4})$ with $\sigma=0.02$, $\psi^{(1)}(x,y)=0$
\item Boundary conditions are Dirichlet conditions consistent with the initial conditions. Because the signal does not reach the boundary, these conditions do not influence the result.
\end{itemize}

Due to the rotational invariance of the initial condition and the flux of the $M_1$ model the exact solution for the line-source problem has rotational symmetry.

\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{Linesource4}{\input{Images/Linesource4}}
\caption{Line source: comparison of primitive and characteristic slope limiter.}
\label{fig:linesource4}
\end{figure}

This symmetry is easily destroyed by the application of the slope limiter in primitive variables (SL$0$$\square/\triangle$, Figure~\ref{fig:linesource4}). 
It is visible for both mesh types that the solution is strongly asymmetric in density as well as normalized velocity. Furthermore, realizability is lost in parts of the domain (black and white spots).

Applying the characteristic limiter (CL$0$$\square/\triangle$) yields much better results: the solution is much more symmetric and as shown in 
Tables~\ref{tab:linesource_rect_table} and~\ref{tab:linesource_tri_table}, it is even realizable in the mean at the final time step. 
However, the realizability in the mean does not hold at every time step and in addition there are a Gauss points on both meshes at which the solution is unrealizable. 
We note that the fact that primitive-variable limiting (SL$0$$\square/\triangle$) can result in non-physical oscillations is a well-known phenomenon 
and has been reported in case of the $M_1$ model in one dimension in \cite{Olbrant12}. 

\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{Linesource4b}{\input{Images/Linesource4b}}
\caption{Line source: comparison of the characteristic limiter with and without realizability limiter.}
\label{fig:linesource4b}
\end{figure}

Unfortunately, choosing $M=0$ in the minmod limiter results in a flattening of the solution at smooth extrema and first-order convergence in the $L_\infty$-norm \cite{CockburnShuPk}. 
Setting $M=22$ results in sharper solutions in the sense that the width of the wave front is decreasing (Figure \ref{fig:linesource4b}). However, this comes at the price of a realizability loss resulting in a complete destruction of the solution (see e.g. CL$22\square$ in Figure~\ref{fig:linesource4b}). Similar results occur with the primitive-variable limiter (not shown).
This can be fixed by further application of the realizability limiter (CRL$22\square/\triangle$). The now-realizable solutions show the same sharpness as with CL$22\square/\triangle$. 
Furthermore, good rotational symmetry is still visible.

\begin{figure}[!h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{Linesource2}{\input{Images/Linesource2}}
\caption{Line source: realizability limiter only.}
\label{fig:linesource2}
\end{figure}

A natural question arising is if the realizability limiter on its own is enough to provide a good solution. It has been shown in \cite{Olbrant12} that even in one dimension spurious oscillations can occur. Similarly, as shown in Figure \ref{fig:linesource2}, symmetry is not preserved and slight oscillations are visible, especially in the normalized velocity. Although the solution is everywhere realizable, a slope limiter is strictly necessary to dampen the Gibbs phenomenon. Similar observations have been made in case of the Euler equations in \cite{Zhang2010,Wang2012}.

\begin{figure}[h]
\def\localpath{./Images/Linesource/}
\centering
\externaltikz{Linesource4c}{\input{Images/Linesource4c}}
\caption{Line source: primitive variable and realizability limiter.}
\label{fig:linesource4c}
\end{figure}

As we have shown, using the characteristic and realizability limiter together yields good results. However, the change to characteristic fields is expensive and the transformation matrices become arbitrarily poorly conditioned close to the realizability boundary \cite{AllHau12}, yielding inaccurate results. Unfortunately, Figure \ref{fig:linesource4c} shows that the combination of limiting in the 
primitive-variable and realizability limiter is not enough to obtain good results. It appears that these solutions are even worse than the realizability-limiter-only solutions shown in Figure~\ref{fig:linesource2}.

    Our numerical tests show that in general the characteristic limiter combined with the realizability limiter (CRLM) is superior
    to the combination of limiting in the components and the realizability limiter (SRLM). However, the differences between the two limiting
    strategies are most dramatic for the linesouce test case. In the rest of the problems we consider, the solutions both look good
    and are comparable.

\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\square$}& \multicolumn{2}{c}{SL$0\square$}& \multicolumn{2}{c}{SL$22\square$} & \multicolumn{2}{c}{CL$0\square$}& \multicolumn{2}{c}{CL$22\square$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP & 11 & 49\%& 8 & 19\%& 13 & 01\%& 0 & 84\%& 7 & 69\%\\
$\max\limits_t$ CM & 10 & 09\%& 8 & 08\%& 10 & 59\%& 0 & 81\%& 5 & 87\%\\
GP$(t_{f})$ & 11 & 49\%& 3 & 26\%& 13 & 01\%& 0 & 06\%& 7 & 69\%\\
CM$(t_{f})$ & 10 & 09\%& 3 & 15\%& 10 & 59\%& 0 & 00\%& 5 & 87\%\\
\end{tabular}
\caption{Line source: percentage of non-realizable cells and quadrature points on rectangular mesh.}
\label{tab:linesource_rect_table}
\end{table}
\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\triangle$}& \multicolumn{2}{c}{SL$0\triangle$}& \multicolumn{2}{c}{SL$22\triangle$} & \multicolumn{2}{c}{CL$0\triangle$}& \multicolumn{2}{c}{CL$22\triangle$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP & 81 & 01\%& 5 & 01\%& 30 & 10\%& 0 & 03\%& 0 & 33\%\\
$\max\limits_t$ CM & 80 & 33\% & 5 & 03\%& 25 & 12\%& 0 & 00\%& 0 & 17\%\\
GP$(t_{f})$ & 81 & 01\%& 5 & 05\%& 30 & 10\% & 0& 01\%& 0& 20\%\\
CM$(t_{f})$ & 80 & 33\% & 5 & 00\%& 25 & 12\%& 0 & 00\%& 0 & 08\%\\
\end{tabular}
\caption{Line source: percentage of non-realizable cells and quadrature points on triangular mesh.}
\label{tab:linesource_tri_table}
\end{table}

In Tables~\ref{tab:linesource_rect_table} and \ref{tab:linesource_tri_table} we summarize the performance of the slope limiter in the absence of the realizability 
limiter by counting the proportion of non-realizable cells (CM) and Gauss points (GP). We provide both the maximum values in the time interval (taken over a sample of $20$ equidistant
time points) and the values at the final time step for each mesh. It is clear that the characteristic limiter (CL$M$$\square/\triangle$) performs better than the slope limiter in primitive variables
(SL$M$$\square/\triangle$). In fact the final solution from the characteristic limiter (CL$0$$\square/\triangle$) is realizable in the mean on both meshes. However, realizability in the mean and at the Gauss points is not guaranteed
at all time steps on the time interval. These results further reinforce the need for the realizability limiter.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Hdisk
%\newpage
\subsection{Homogeneous Disk}
Our next test case is the two-dimensional version of the homogeneous sphere test, which is often used to test radiation transport codes, see e.g. \cite{Abdikamalov2012,Bruenn1985,Rampp2002,Sumiyoshi2012}. It consists of a static homogeneous and isothermal sphere that radiates in vacuum \cite{Radice2013648}. It is well-suited to test moment approximations since it admits an analytical solution \cite{smit1997hyperbolicity}. The parameters in the two-dimensional setting are as follows:
\begin{itemize}
\item Domain: $[-5,5]^2$ with $h=0.05$
\item Time: 3.00
\item Parameters: Let $D$ be the unit disk $D=\{ (x,y):\ x^2+y^2\leq 1\}$. Set $\sigma_s=0$; $\sigma_a=10$ on $D$, zero otherwise; source $q^{(0)}=1$ on $D$, zero otherwise; $q^{(1)} = 0$
\item Initial condition: $\psi^{(0)}=10^{-10}$, $\psi^{(1)}=0$
\item Boundary conditions: $\psi^{(0)}=10^{-10}$, $\psi^{(1)}=0$
\end{itemize}

This test case is numerically challenging due to the discontinuous parameters. Furthermore, its solution is again radially symmetric, simplifying the analysis of the quality of the limiter configurations.

\begin{figure}[h]
\def\localpath{./Images/HomDisk/}
\centering
\externaltikz{Hdisk4a}{\input{Images/Hdisk4a}}
\caption{Homogeneous disk: comparison of primitive and characteristic slope limiter.}
\label{fig:Hdiska}
\end{figure}
As in the line source test case, we observe that the application of the primitive-variable limiter results in an oscillatory solution which is not realizable. This is shown in Figure~\ref{fig:Hdiska}. In case of the triangular mesh strong filaments form at the edge of the disk. Using the characteristic limiter slightly improves the solution. However, and in contrast to the line source, even with $M=0$ 
the characteristic limiter is no longer able to preserve realizability in the mean (see Tables~\ref{tab:hdisk_rect_table} and~\ref{tab:hdisk_tri_table}), furthermore, a larger proportion of
the Gauss points are not realizable. Additionally, at the boundary of the disk the quality of the solution is strongly reduced due to the bad condition number of the transformation matrix, which 
goes to infinity for values of $\lvert\phi^{1}\rvert$ approaching $1$ \cite{AllHau12}.

\begin{figure}[h]
\def\localpath{./Images/HomDisk/}
\centering
\externaltikz{Hdisk4b}{\input{Images/Hdisk4b}}
\caption{Homogeneous disk: comparison of the characteristic limiter with and without realizability limiter}
\label{fig:Hdiskb}
\end{figure}
%
Adding the realizability limiter strongly improves the quality of the solution (Figure~\ref{fig:Hdiskb}). As before, choosing $M=0.2$ gives less diffusive results.

\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\square$}& \multicolumn{2}{c}{SL$0\square$}& \multicolumn{2}{c}{SL$0.2\square$} & \multicolumn{2}{c}{CL$0\square$}& \multicolumn{2}{c}{CL$0.2\square$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP  & 24 & 13\%& 12 & 46\%& 18 & 95\%& 1 & 57\%& 5 & 53\%\\
$\max\limits_t$ CM  & 21 & 95\%& 7 & 19\%& 15 & 08\%& 0 & 22\%& 1 & 14\%\\
GP$(t_{f})$ & 4& 93\%&12& 46\%& 18 & 46\%& 1 & 23\%& 2 & 91\%\\
CM$(t_{f})$ & 0& 26\%& 7 & 19\%& 9 & 39\%& 0 & 06\%& 0 & 06\%\\
\end{tabular}
\caption{Homogeneous disk: percentage of non-realizable cells and quadrature points on rectangular mesh.}
\label{tab:hdisk_rect_table}
\end{table}
%
\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\triangle$}& \multicolumn{2}{c}{SL$0\triangle$}& \multicolumn{2}{c}{SL$0.2\triangle$} & \multicolumn{2}{c}{CL$0\triangle$}& \multicolumn{2}{c}{CL$0.2\triangle$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP & 71 & 75\%& 34&25\% &  41 & 73\%& 5 & 65\% & 33 & 65\%\\
$\max\limits_t$ CM & 71 & 62\% & 34&41\% & 40 & 18\%& 5 & 75\% & 32 & 64\%\\
GP$(t_{f})$ & 71 & 75\%& 34 & 25\%& 41 & 73\%& 5 & 65\%& 33 & 65\%\\
CM$(t_{f})$ & 71 & 61\%& 34 & 41\%& 40 & 18\%& 5 & 71\%& 32 & 64\%\\
\end{tabular}
\caption{Homogeneous disk: percentage of non-realizable cells and quadrature points on triangular mesh.}
\label{tab:hdisk_tri_table}
\end{table}
%
A quantitative analysis of the slope limiters (SL$M$$\square/\triangle$) and (CL$M$$\square/\triangle$) show that limiting in characteristic variables
yields better results but the solution remains unrealizable.
%\newpage
\subsection{Flash}
%\comment{M=0.5}

In this test case a bulk of mass is moving from the center of the domain to the right boundary \cite{Science2013}. The parameters are given as follows:
\begin{itemize}
\item Domain: $[-10,10]^2$ with $h =0.06$
\item Time: 6
\item Parameters: Let $D$ be unit disk with radius $0.5$, $D=\{ (x,y):\ x^2+y^2\leq 0.25\}$. Set $\sigma_s=0$; $\sigma_a=0$; source $q^{(0)}=0$, $q^{(1)} = 0$
\item Initial condition: $\psi^{(0)}=1$, $\psi^{(1)}_x=0.9$, $\psi^{(1)}_y=0$ on $D$, $\psi^{(0)}=10^{-10}$, $\psi^{(1)}=0$ otherwise
\item Boundary conditions: $\psi^{(0)}=10^{-10}$, $\psi^{(1)}=0$
\end{itemize}

Due to the initial condition which is placed close to the boundary of realizability this test is well suited to show the inability of standard slope limiters to preserve realizability. This is demonstrated in Figure~\ref{fig:flash4a} where similar effects as before occur.

\begin{figure}[htbp]
\def\localpath{./Images/Flash/}
\centering
\externaltikz{flash4a}{\input{Images/flash4a}}
\caption{Flash test: comparison of primitive and characteristic slope limiter.}
\label{fig:flash4a}
\end{figure}

\begin{figure}[htbp]
\def\localpath{./Images/Flash/}
\centering
\externaltikz{flash4}{\input{Images/flash4}}
\caption{Flash test: characteristic limiter with realizability limiter for different values of $M$.}
\label{fig:flash4}
\end{figure}

A comparison of Figures~\ref{fig:flash4a} and~\ref{fig:flash4} confirms that the combination of the characteristic limiter with the realizability
limiter yields the best solutions. The performance of the slope limiter without the realizability limiter in Tables~\ref{tab:flash_rect_table} and~\ref{tab:flash_tri_table}
is consistent with previous observations in the line source and homogeneous disk test cases.

\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\square$}& \multicolumn{2}{c}{SL$0\square$}& \multicolumn{2}{c}{SL$0.5\square$} & \multicolumn{2}{c}{CL$0\square$}& \multicolumn{2}{c}{CL$0.5\square$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP & 14 & 84\%& 13 & 97\%& 20 & 88\%& 0 & 86\%& 19 & 75\%\\
$\max\limits_t$ CM & 11 & 77\%& 11 & 53\%& 14 & 27\%& 0 & 22\%& 15 & 01\%\\
GP$(t_{f})$ & 14 & 70\%& 13 & 97\%& 20 & 88\%& 0 & 86\%& 19 & 75\%\\
CM$(t_{f})$ & 9 & 50\%& 11 & 53\%& 14 & 27\%& 0 & 22\%& 15 & 01\%\\
\end{tabular}
\caption{Flash: percentage of non-realizable cells and quadrature points on rectangular mesh}
\label{tab:flash_rect_table}
\end{table}

\begin{table}
\centering
\begin{tabular}{l|| r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l }
& \multicolumn{2}{c}{SL$\infty\triangle$}& \multicolumn{2}{c}{SL$0\triangle$}& \multicolumn{2}{c}{SL$0.5\triangle$} & \multicolumn{2}{c}{CL$0\triangle$}& \multicolumn{2}{c}{CL$0.5\triangle$}\\
\cmidrule(r){1-1}\cmidrule(r){2-3} \cmidrule(r){4-5} \cmidrule(r){6-7} \cmidrule(r){8-9} \cmidrule(r){10-11} \
$\max\limits_t$ GP & 55 & 21\%& 24 & 12\%& 50 & 14\%& 5 & 45\%& 50 & 99\%\\
$\max\limits_t$ CM & 55 & 18\% & 24 & 27\% & 49 & 10\% & 5 & 63\%& 39 & 00\%\\
GP$(t_{f})$ & 55 & 28\%& 24 & 12\%& 50 & 14\%& 5 & 45\%& 40 & 99\%\\
CM$(t_{f})$ & 55 & 18\%& 24 & 27\%& 49 & 10\%& 5 & 63\%& 39 & 00\%\\
\end{tabular}
\caption{Flash: percentage of non-realizable cells and quadrature points on triangular mesh.}
\label{tab:flash_tri_table}
\end{table}
outlook.tex
\section{Conclusions and outlook}
In this work, we have investigated a third order realizability-preserving DG method to approximate solutions to the $M_1$ model of radiation transport. The results show that in all test cases presented, the realizability limiter is imperative to ensure the realizability of the numerical solutions.

Based on these results, we advocate the $M_1$ model as an ideal test case that pushes the boundaries of high-order realizability-preserving numerical methods. The model equations are highly nonlinear, and are well-posed only in a geometrically complex (yet convex) domain of realizability. Furthermore, practically relevant solutions are not just wave fronts, but rather beam-like solutions, or they possess Dirac-like source terms. Solutions and material coefficients tend to vary over several orders of magnitude. Furthermore, so-called source-detector problems in radiation transport require an accurate solution (in terms of the relative error) at one or more specific points in space.

The results presented in this paper cast some doubt on whether DG is a competitive numerical scheme for $M_1$, considering the accuracy requirements commonly encountered in practice. For instance, for the homogeneous disk test case the DG solution possesses around $10^6$ degrees of freedom. It is comparable in quality to a Lax-Friedrichs (LF) solution with roughly the same number of degrees of freedom ($2048\times 2048\approx 4\cdot 10^6$), which however is significantly faster to compute because no limiting is needed. Similarly, in the line source test case, the best DG solution is comparable to a $1024\times 1024$ LF solution.

In several test cases, the high-order DG method yields reliable results only if a full limiting strategy is employed, including transformations to and from the characteristic variables. Not only do those transformations incur a significant cost, they also become less and less well conditioned when higher accuracy requirements are imposed. It is therefore questionable whether, for the $M_1$ model, the high-order DG methodology will ever become more efficient than a simple monotone first-order scheme that automatically guarantees realizability.

That being said, it should be pointed out that the DG methodology studied here, albeit being an established standard method, can most certainly be improved. For instance, more sophisticated limiting techniques, such as WENO limiting with KXRCF shock detection \cite{Krivodonova2004,Qiu2005a}, or limiters that incorporate high-order approximations during the limiting process~\cite{Burbeau2001,Krivodonova2007}, could be combined with the realizability-preserving limiter. This would remove the need for the problem-dependent minmod parameter $M$. Furthermore, more efficient time stepping schemes could be applied, such as explicit low-storage SSP Runge-Kutta schemes \cite{Ketcheson2008}, SSP multi-step Runge-Kutta methods \cite{Ketcheson,Bresten2013} (which both require less evaluations of the differential operator $\mathcal{L}_h$), or even Implicit-Explicit schemes \cite{Ascher1997,Kanevsky2007} (to remove the dependence on the right-hand side ($\sigma_a$ and $\sigma_s$) in the CFL condition \eqref{eqn:cfl}).

Further research steps include realizability limiting for higher-order minimum-entropy models (i.e., $M_N$ with $N\geq 2$) in two and three dimensions, similar to \cite{Schneider2015a}. The main challenge in that context is that a complete characterization of the realizability set by means of simple algebraic criteria is still outstanding. One possible remedy is to use discrete realizability criteria with respect to a quadrature rule \cite{AllHau12}. Moreover, in general a high number of quadrature points is necessary to produce solutions without discernible discrete ray artifacts \cite{Schneider2016,Schneider2015a,Garrett2014}, resulting in an extremely high cost for realizability limiting. Hence less expensive strategies to guarantee realizability need to be developed.




properties.tex
\section{Properties of the $M_1$ Model}
\label{sec:M1}
%
From the general theory of moments one can deduce
conditions on a given set of values, that are necessary and sufficient
for the existence of a probability measure whose moments match these
values \cite{Ker76}. Such moments are called realizable. In our case, we can restrict ourselves to measures that have a formal density (i.e.\ for the closure we allow Dirac $\delta$'s as densities).
Given a scalar $\psiz$ and a vector $\psio$, these are the zeroth and first moment of a non-negative density if and only if \cite{Ker76,Schneider2016,Monreal}
$$
\psiz>0 \text{ and } |\psio|\leq \psiz.
$$
The interior of the realizable set is described by $\psiz>0 \text{ and } |\psio| < \psiz$.
It has been shown in \cite{Goudon2005} that, under reasonable assumptions on the initial conditions, the analytical solution of the $M_1$ model in one spatial dimension remains realizable for all time. Until now, no similar statement could be made for the $M_1$ model in multiple dimensions or higher-order $M_N$ models. 

As mentioned before, the hyperbolicity of the $M_1$ model is a direct consequence of the entropy closure. Additional insight can be obtained by computing the eigenvalues explicitly. Given a unit direction vector $\bfn=(n_x,n_y)^T$, we have to determine the eigenvalues of the directional Jacobian
\begin{equation}
\label{eq:jacobian}
n_xJ_x(U) + n_yJ_y(U),
\end{equation}
where $J_x$ and $J_y$ are the Jacobians of $F$ and $G$, respectively. Because of the invariance of the model under coordinate transformations, the eigenvalues can only depend on the angle between $\psio$ and $\bfn$, and the absolute value $\frac{|\psio|}{\psiz}$. The eigenvalues can be computed analytically, but the formulas are very lengthy so we do not show them here (cf.\ \cite{Berthon2006}). Figure \ref{fig:EV} shows the eigenvalues for two different angles. The most important observation is that the eigenvalues collapse into one value at the boundary of the realizability region, i.e.\ for $\frac{|\psio|}{\psiz}=1$. Inspection of the Jacobian shows that at that boundary the Jacobian is no longer diagonalizable, i.e.\ the $M_1$ model loses (strict) hyperbolicity. This means that the $M_1$ model is only well-posed in the interior of the domain of realizability. As a consequence,  in a numerical scheme one should always ensure staying in the interior of the relizability domain.

\begin{figure}

\centering
\externaltikz{eigenvalues1}{\input{Images/M1eigenvalues1}}  
\caption{Eigenvalues of the $M_1$ model as a function of $f=\frac{n\cdot\psio}{\psiz}$ for a fixed angle $\varphi$ between $\psio$ and $n$.}
\label{fig:EV}
\end{figure}
realizabilitylimiter.tex
\section{Realizability Limiter}
\label{sec:limiter}
We want to construct a scheme so that the numerical solution stays realizable with respect to the cell means. The proofs presented in this section follow the strategy used in {\cite{Zhang2010}, which has also been used in the construction of realizability-preserving limiters for the $M_1$ model in 1D \cite{Olbrant12}. 

First we need a technical lemma.
\begin{lem}
\label{lem:Lemma3}
Let $\nu\in \mathbb{R}^2$ be an arbitrary unit vector. For the $M_1$ model the combination of moments $(\psiz\pm \psio\cdot\nu,\psio\pm \psit\cdot \nu)$ is realizable.
\end{lem}
\begin{proof}
Let $\psi$ be a non-negative density that realizes $\psiz$ and $\psio$. Then $(1\pm\nu\cdot\Omega)\psi(\Omega)$ is a non-negative density that realizes the combined moments.
\end{proof}

Let us now consider the higher-order scheme \eqref{eq:dweakform2}. Due to the convexity of the realizable set it suffices to investigate the forward-Euler discretization in time since the used SSP integrator is just a convex combination of such Euler steps.
\begin{thm}[Main Result]
\label{thm:MainResult}
One forward-Euler step of the scheme \eqref{eq:dweakform2} with the DG polynomial $p_K$ of degree $k$ yields realizable cell averages
if $p_K(x_s,y_s)$ is realizable for all $(x_s,y_s)\in S_k^K$ and if the CFL condition
\begin{gather}
\frac{2}{3}\hat{w}_1\left(1-\Delta t\left(\sigma_a + \sigma_s\right)\right) - \frac{\Delta t}{2\abs{K}}l_K^i\geq 0~~\forall i=1,\ldots,3
\label{eqn:cfl}
\end{gather}
holds. Here $\hat{w}_1$ is the quadrature weight of the $N$-point Gauss-Lobatto rule on $[-\frac{1}{2},\frac{1}{2}]$ for the first quadrature point.
\end{thm}
\begin{proof}
    For simplicity assume that $\sigma_s=\sigma_a = 0$ and $q = 0$. Furthermore, to ease notation, we drop the symbol $n$ for all quantities at the current time step. After application of the forward-Euler scheme the cell averages in element $K$ satisfy 
\begin{align*}
\overline{U}_K^{n+1} &= \overline{U}_K - \frac{\Delta t}{\abs{K}}\sum\limits_{i=1}^{3}\int_{e_i} H_{e_i,K}\left(U_{K}^{int},U_{K(i)}^{ext}\right)~dx\\
&\stackrel{\eqref{eq:LFFlux}}{=}  \begin{pmatrix}\psizbar_K\\\psiobar_K\end{pmatrix} - \frac{\Delta t}{2\abs{K}}\sum\limits_{i=1}^{3}\int_{e_i} \begin{pmatrix}\psiz_K\\\psio_K\end{pmatrix}~dx + \frac{\Delta t}{2\abs{K}}\sum\limits_{i=1}^{3}\int_{e_i} \begin{pmatrix}\psiz_{K(i)}-\psio_{K(i)}\cdot \bfn_{e_i,K}\\\psio_{K(i)} - \psit_{K(i)}\cdot \bfn_{e_i,K}\end{pmatrix}~dx,
\end{align*}
where $U_K = p_K$ and $U_{K(i)} = p_{K(i)}$ denote the polynomial representations of $U_h$ on element $K$ and on its neighbor $K(i)$ sharing edge $e_i$ with $K$, respectively.

Evaluating the second integral using the Gaussian quadrature as in \eqref{eq:FluxGauss} it is easy to see by using Lemma \ref{lem:Lemma3} that on each edge $e_i$ the quantity
$$
\int_{e_i} \begin{pmatrix}\psiz_{K(i)}-\psio_{K(i)}\cdot \bfn_{e_i,K}\\\psio_{K(i)} - \psit_{K(i)}\cdot \bfn_{e_i,K}\end{pmatrix}~dx 
=  
l_K^i\sum_{\beta=1}^{k+1} w^\beta \begin{pmatrix}\psiz_{K(i),\beta}-\psio_{K(i),\beta}\cdot \bfn_{e_i,K}\\\psio_{K(i),\beta} - \psit_{K(i),\beta}\cdot \bfn_{e_i,K}\end{pmatrix}
$$
is realizable since $\sum w^\beta = 1$. Here, $\psiz_{K(i),\beta}$ denotes the evaluation of the polynomial representation of $\psiz$ on the neighboring element $K(i)$ at the quadrature node associated with $w^\beta$ on $e_i$ and similarly for all other moments.

Now it suffices to show that the first terms in this equation are realizable, too. We rewrite the cell average into a combination of interior points and edge points:
\begin{align*}
\overline{U}_K &\stackrel{\eqref{eq:Cellmean}}{=} \sum\limits_{i=1}^{3}\sum\limits_{\beta=1}^{k+1}\frac{2}{3}w^\beta\hat{w}^1 U_{i,\beta}^{int} + \sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma U_\gamma^{inner},\\
\int_{e_i} U_K~dx  &= \sum\limits_{\beta=1}^{k+1}w^\beta l_K^i U_{i,\beta}^{int}.
\end{align*}
Then
\begin{align*}
\begin{pmatrix}\psizbar_K\\\psiobar_K\end{pmatrix} - \frac{\Delta t}{2\abs{K}}\sum\limits_{i=1}^{3}\int_{e_i} \begin{pmatrix}\psiz_K\\\psio_K\end{pmatrix}~dx = \sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma U_\gamma^{inner} + \sum\limits_{i=1}^{3}\sum\limits_{\beta=1}^{k+1}w^\beta U_{i,\beta}^{int} \left(\frac{2}{3}\hat{w}^1 - \frac{\Delta t}{2\abs{K}}l_K^i\right).
\end{align*}
Under the CFL condition, we immediately get that 
$$
\sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma U_\gamma^{inner} + \sum\limits_{i=1}^{3}\sum\limits_{\beta=1}^{k+1}w^\beta U_{i,\beta}^{int} \left(\frac{2}{3}\hat{w}^1 - \frac{\Delta t}{2\abs{K}}l_K^i\right)
$$
is realizable since the quadrature weights sum up to $1$.\\ 
For absorption and scattering we assume for simplicity that $\sigma_s$ and $\sigma_a$ are constant in $K$. Then we can write the cell averages of the source term as 
\begin{align*}
\overline{S(U_K)} &=  \sum\limits_{i=1}^{3}\sum\limits_{\beta=1}^{k+1}\frac{2}{3}w^\beta\hat{w}^1 \left(\sigma_a U_{i,\beta}^{int} + \sigma_s \begin{pmatrix}
0\\\psio_{i,\beta,int}
\end{pmatrix} \right)\\&+ \sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma \left(\sigma_a U_\gamma^{inner} + \sigma_s \begin{pmatrix}
0\\\psio_{\gamma,{inner}}
\end{pmatrix} \right).
\end{align*}
Thus the updated cell averages can be written as the moments of a non-negative distribution, i.e. $U_K^{n+1} = \int \begin{pmatrix}
1\\\Omega
\end{pmatrix}\Xi_K(\Omega)~d\Omega$ with
\begin{align*}
\Xi_K &=\sum\limits_{\gamma=1}^{L}\tilde{w}^\gamma \left[\left(1-\Delta t\left(\sigma_a + \sigma_s\right)\right)\psi_\gamma^{inner} + \cfrac{\sigma_s}{4\pi}\int \psi_\gamma^{inner}~d\Omega \right] \\
&+ \sum\limits_{i=1}^{3}\sum\limits_{\beta=1}^{k+1}w^\beta \left[\psi_{i,\beta}^{int} \left(\frac{2}{3}\hat{w}^1\left(1-\Delta t\left(\sigma_a + \sigma_s\right)\right) - \frac{\Delta t}{2\abs{K}}l_K^i\right) +\frac{\Delta t}{2\abs{K}}l_K^i {\psi}_{K(i),\beta}+ \cfrac{\sigma_s}{4\pi}\int \psi_{i,\beta}^{int} ~d\Omega\right]
\end{align*}
The quantity $\psi_\gamma^{inner}$ is any nonnegative distribution function realizing $U_{\gamma}^{inner}$, similarly for $\psi_{i,\beta}^{int}$ and ${\psi}_{K(i),\beta}$. It is easy to see that under the given assumptions $\Xi_K^n\geq 0$, which implies that $U_K^{n+1}$ is by definition realizable.
\end{proof}

A similar result can be obtained on the rectangular grid by going through the same lines of the previous proof, replacing the quadrature rules on the triangle with those on the rectangles. More details can be found in \cite{Schneider2016}.

\begin{thm}[Main Result for rectangular grid]
\label{thm:MainResultRect}
One forward-Euler step of the scheme \eqref{eq:dweakform2} with the DG polynomial $p_K$ of degree $k$ yields realizable cell averages
if $p_K(x_s,y_s)$ is realizable for all $(x_s,y_s)\in S_k^K$ and if the CFL condition
$$
\hat{w}_1\left(1-\Delta t\left(\sigma_a + \sigma_s\right)\right) - \frac{\Delta t}{\Delta x}- \frac{\Delta t}{\Delta y}\geq 0
$$
holds. Here $\hat{w}_1$ is the quadrature weight of the $N$-point Gauss-Lobatto rule on $[-\frac{1}{2},\frac{1}{2}]$ for the first quadrature point.
\end{thm}

All that remains is to ensure that the assumptions of the two theorems are satisfied. This can be achieved by a simple scaling limiter. A similar limiter has been derived in other contexts, e.g. shallow water \cite{Xing2013}, Euler equations \cite{Zhang2010} and gas dynamics \cite{Zhang2012a}. The basic idea is to dampen the higher-order parts of the DG polynomial $p_K(x,y)$ until at all quadrature nodes $(x_s,y_s)\in S_k^K$ the evaluation $p_K(x_s,y_s)$ is realizable. It has been reported in \cite{Zhang2012a} that this limiter, due to its simplicity, can destroy the formal accuracy of the scheme in certain non-generic situations. This has been investigated in a one-dimensional setting for minimum-entropy models in \cite{Schneider2015a}. There, it has been shown that close to the realizability boundary the convergence order can drop.  However, we still expect order preservation in most practical cases.

Writing the limited polynomial as 
\begin{align*}
p_K^\theta(x,y) := \theta \overline{p_K} + (1-\theta)p_K(x,y) \stackrel{\eqref{eqn:solution_form}}{=} \theta \alpha_1^U\varphi_1 +  (1-\theta)\sum\limits_{i=1}^{N_k} \alpha_i^U\varphi_i = \alpha_1^U\varphi_1 + (1-\theta)\sum\limits_{i=2}^{N_k} \alpha_i^U\varphi_i
\end{align*}
it becomes clear that the cell average of $p_K$ is preserved and only higher-order coefficients are damped by the factor $(1-\theta)$ where $\theta$ is chosen as the minimal value in the set
\begin{align}
\left\{\tilde{\theta}\in[0,1]~|~p_K^{\tilde\theta}(x_s,y_s)~\text{ is realizable for all } (x_s,y_s)\in S_k^K \right\}.
\label{eq:RLproblem}
\end{align}
Such a $\theta$ always exist under the assumption that $\overline{p_K}$ is realizable as a consequence of the convexity of the realizable set. To calculate $\theta_s$ it suffices to find the intersections of the ray segments
\begin{align}
\label{eq:RaySegment}
\left\{\theta \bar U_K +(1-\theta) p_K(x_s,y_s)~|~\theta\in[0,1]\right\},~(x_s,y_s)\in S_k^K
\end{align} 
with the boundary of the realizable set
\begin{align}
\label{eq:BoundaryReal}
\left\{\left(\psi^{(0)},\psi^{(1)}\right)\in\mathbb{R}^3~|~ \psi^{(0)} = \left\lvert \psi^{(1)}\right\rvert \right\} = \left\{\left(\psi^{(0)},\psi^{(1)}\right)\in\mathbb{R}^3~|~ \psi^{(0)}\geq 0, \left(\psi^{(0)}\right)^2 = \left\lvert\psi^{(1)}\right\rvert^2 \right\}.
\end{align}
Using the condition $\left(\psi^{(0)}\right)^2 = \left\lvert\psi^{(1)}\right\rvert^2$ from \eqref{eq:BoundaryReal} and the parameterization of the ray segment \eqref{eq:RaySegment} we get
\begin{align*}
\left(\theta\bar\psi^{(0)}+(1-\theta)\psi^{(0)}\right)^2 \stackrel{!}{=} \left(\theta\bar\psi^{(1)}_y+(1-\theta)\psi^{(1)}\right)^2+\left(\theta\bar\psi^{(1)}_y+(1-\theta)\psi^{(1)}\right)^2.
\end{align*} 
Rearranging this expression, collecting the coefficients in powers of $\theta$ yields the polynomial equation
\begin{align*}
0\stackrel{!}{=} r(\theta) &= a\theta^2+b\theta+c,\\
a &= -\left(\bar{\psi}^{(0)}-{\psi}^{(0)}\right)^2 + \left(\bar{\psi}^{(1)}_x -{\psi}^{(1)}_x\right)^2 + \left(\bar{\psi}^{(1)}_y -{\psi}^{(1)}_y\right)^2\\
b &= 2\cdot\left(\left({\psi}^{(0)}\right)^2 - \bar{\psi}^{(0)}\cdot{\psi}^{(0)} - \left({\psi}^{(1)}_x\right)^2 + \bar{\psi}^{(1)}_x\cdot{\psi}^{(1)}_x - \left({\psi}^{(1)}_y\right)^2 + \bar{\psi}^{(1)}_y\cdot{\psi}^{(1)}_y\right)\\
c &= - \left({\psi}^{(0)}\right)^2 + \left({\psi}^{(1)}_x\right)^2 + \left({\psi}^{(1)}_y\right)^2.
\end{align*}
Since $r(\theta)$ is a quadratic polynomial, it possibly has two zeros in the set $[0,1]$. Due to convexity, the largest zero in this set is the correct choice for $\theta_s$, since it corresponds to the intersection with the $\psi^{(0)}\geq 0$ part of the double cone $|\psi^{(1)}|^2\leq (\psi^{(0)})^2$, while the smaller zero (which is by definition farther away from $\bar U_K$) corresponds to the intersection with the  $\psi^{(0)}\leq 0$ part of the double cone. Finally, set $\theta := \max_s \theta_s$.

Note that in each intermediate step of the Runge-Kutta scheme the realizability limiter is always applied after the slope limiting. Applying the slope limiter can destroy the realizability again \cite{Schneider2016}.








usermacros.tex
\newcommand{\prn}[1]{\left(#1\right)}
\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\brk}[1]{\left<#1\right>}
\newcommand{\pair}[2]{\left<#1,#2\right>_{SS'}}
\newcommand{\cR}{\mathcal{R}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\mean}[1]{\langle #1 \rangle}
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\mat}[1]{\mathbf{#1}}

\newcommand{\eps}{\varepsilon}
\newcommand{\Psihat}{\widehat\Psi}
\newcommand{\psihat}{\widehat\psi}
\newcommand{\phihat}{\widehat\phi}
\newcommand{\rhohat}{\widehat\rho}
\newcommand{\Qhat}{\widehat Q}

\newcommand{\bfU} { \mbox{\boldmath $U$} }
\newcommand{\bfuz} { \mbox{\boldmath $u_0$} }
\newcommand{\bfG} { \mbox{\boldmath $G$} }
\newcommand{\bfn} { \mbox{\boldmath $n$} }
\newcommand{\bfF} { \mbox{\boldmath $F$} }
\newcommand{\bfM} { \mbox{\boldmath $M$} }
\newcommand{\bfH} { \mbox{\boldmath $H$} }
\newcommand{\bfS} { \mbox{\boldmath $S$} }
\newcommand{\sF}  { \mbox{$\mathcal{F}$} }
\newcommand{\bfalpha} { \mbox{\boldmath $\alpha$} }
\newcommand{\bfgamma} { \mbox{\boldmath $\gamma$} }


\newcommand{\psiz}{\psi^{(0)}}
\newcommand{\psio}{\psi^{(1)}}
\newcommand{\psit}{\psi^{(2)}}
\newcommand{\psizbar}{\bar{\psi}^{(0)}}
\newcommand{\psiobar}{\bar{\psi}^{(1)}}

\pgfplotsset{
plotstylea/.style={royalblue!70,every mark/.append style={solid,line width = 0pt,fill=royalblue!60!black},mark=ball}, plotstyleb/.style={dashed,no marks}, }

DGforM1.bbl
\begin{thebibliography}{10}

\bibitem{Abdikamalov2012}
E.~Abdikamalov, A.~Burrows, C.~D. Ott, F.~L{\"{o}}ffler, E.~O'Connor, J.~C.
  Dolence, and E.~Schnetter.
\newblock {A New Monte Carlo Method for Time-Dependent Neutrino Radiation
  Transport}.
\newblock \emph{Astrophys. J.}, 755(2):page 111, 2012.

\bibitem{AllHau12}
G.~Alldredge, C.~Hauck, and A.~Tits.
\newblock High-order entropy-based closures for linear transport in slab
  geometry {II}: A computational study of the optimization problem.
\newblock \emph{SIAM. J. Sci. Comput}, 34(4): pp.  B361--B391, 2012.

\bibitem{Schneider2015a}
G.~W. Alldredge and F.~Schneider.
\newblock {A realizability-preserving discontinuous Galerkin scheme for
  entropy-based moment closures for linear kinetic equations in one space
  dimension}.
\newblock \emph{J. Comput. Phys.}, 295: pp.  665--684, 2015.

\bibitem{Ascher1997}
U.~Ascher, S.~Ruuth, and R.~Spiteri.
\newblock {Implicit--explicit Runge-Kutta methods for time-dependent partial
  differential equations}.
\newblock \emph{Applied Numerical Mathematics}, 1997.

\bibitem{Berthon2006}
C.~Berthon, P.~Charrier, and B.~Dubroca.
\newblock An {HLLC} scheme to solve the ${M}_1$ model of radiative transfer in
  two space dimensions.
\newblock \emph{J. Sci. Comput.}, 31(3): pp.  347--389, 2006.

\bibitem{biswas1994parallel}
R.~Biswas, K.~D. Devine, and J.~E. Flaherty.
\newblock Parallel, adaptive finite element methods for conservation laws.
\newblock \emph{Appl. Numer. Math.}, 14(1): pp.  255--283, 1994.

\bibitem{Bresten2013}
C.~Bresten, S.~Gottlieb, and Z.~Grant.
\newblock {Strong Stability Preserving Multistep Runge-Kutta Methods}.
\newblock \emph{arXiv preprint}, 2013.

\bibitem{Bruenn1985}
S.~W. Bruenn.
\newblock {Stellar core collapse - Numerical model and infall epoch}.
\newblock \emph{Astrophys. J. Suppl. Ser.}, 58:page 771, 1985.

\bibitem{BruHol01}
T.~A. Brunner and J.~P. Holloway.
\newblock One-dimensional {R}iemann solvers and the maximum entropy closure.
\newblock \emph{J. Quant. Spectrosc. Radiat. Transfer}, 69: pp.  543--566,
  2001.

\bibitem{Brunner2005b}
---.
\newblock {Two-dimensional time dependent Riemann solvers for neutron
  transport}.
\newblock \emph{J. Comput. Phys}, 210(September 2001): pp.  386--399, 2005.

\bibitem{Burbeau2001}
A.~Burbeau, P.~Sagaut, and C.-H. Bruneau.
\newblock A problem-independent limiter for high-order runge–kutta
  discontinuous galerkin methods.
\newblock \emph{Journal of Computational Physics}, 169(1): pp.  111 -- 150,
  2001.
\newblock ISSN 0021-9991.

\bibitem{Chidyagwai2016}
P.~Chidyagwai, F.~Schneider, and P.~Monreal.
\newblock {Realizability-preserving Runge-Kutta Discontinuous-Galerkin schemes
  for the M1 model of radiative transfer: Code}, 2016.

\bibitem{Cockburn}
B.~Cockburn.
\newblock {An introduction to the discontinuous Galerkin method for
  convection-dominated problems}.
\newblock \emph{Advanced numerical approximation of nonlinear {\ldots}}.

\bibitem{CockburnShuIV}
B.~Cockburn, S.~Hou, and C.-W. Shu.
\newblock {TVB} {R}unge-{K}utta local projection discontinuous {G}alerkin
  finite element method for conservation laws {IV}: {T}he multi-dimensional
  case.
\newblock \emph{Math. Comp.}, 54: pp.  545--581, 1990.

\bibitem{CockburnLinShu}
B.~Cockburn, S.-Y. Lin, and C.-W. Shu.
\newblock {TVB} {R}unge-{K}utta local projection discontinuous {G}alerkin
  finite element method for consevation laws {III}: {O}ne-dimensional systems.
\newblock \emph{J. Comput. Phys.}, 84: pp.  90 -- 113, 1989.

\bibitem{CockburnShuPk}
B.~Cockburn and C.-W. Shu.
\newblock {TVB} {R}unge-{K}utta local projection discontinuous {G}alerkin
  finite element method for conservation laws {II}: {G}eneral framework.
\newblock \emph{Math. Comp.}, 52(186): pp.  411--435, 1989.

\bibitem{CockburnShuP1}
B.~Cockburn and C.~W. Shu.
\newblock The {R}unge-{K}utta local projection ${P}^1$-discontinuous-{G}alerkin
  finite element method for scalar conservation laws.
\newblock \emph{ESAIM. Math. Model. Num.}, 25(3): pp.  337--361, 1991.

\bibitem{CockburnShuV}
B.~Cockburn and C.-W. Shu.
\newblock The {R}unge--{K}utta discontinuous {G}alerkin method for conservation
  laws {V}: {M}ultidimensional systems.
\newblock \emph{J. Comput. Phys.}, 141(2): pp.  199--224, 1998.

\bibitem{review_article}
---.
\newblock Runge--kutta discontinuous galerkin methods for convection-dominated
  problems.
\newblock \emph{J. Sci.Comput.}, 16(3): pp.  173--261, 2001.

\bibitem{Goudon2005}
J.-F. Coulombel, F.~Golse, and T.~Goudon.
\newblock {Diffusion approximation and entropy-based moment closure for kinetic
  equations}.
\newblock \emph{Asymptotic. Anal}, 45(1,2): pp.  1--39, 2005.

\bibitem{Dre87}
W.~Dreyer.
\newblock Maximisation of the entropy in non-equilibrium.
\newblock \emph{J. Phys. A}, 20: pp.  6505--6517, 1987.

\bibitem{DubFeu99}
B.~Dubroca and J.~L. Feugeas.
\newblock Entropic hierarchy for the radiative transfer equation.
\newblock \emph{C. R. Acad. Sci. Paris Ser. I}, 329: pp.  915--920, 1999.

\bibitem{Frank07}
M.~Frank, H.~Hensel, and A.~Klar.
\newblock A fast and accurate moment method for the {F}okker--{P}lanck equation
  and applications to electron radiotherapy.
\newblock \emph{SIAM. J. Appl. Math.}, 67(2): pp.  582--603, 2007.

\bibitem{Garrett2014}
C.~K. Garrett and C.~D. Hauck.
\newblock {A Comparison of Moment Closures for Linear Kinetic Transport
  Equations: The Line Source Benchmark}.
\newblock \emph{Transp. Theory. Stat. Phys.}, 2013.

\bibitem{gottlieb_shu_tadmor}
S.~Gottlieb, C.-W. Shu, and E.~Tadmor.
\newblock Strong stability-preserving high-order time discretization methods.
\newblock \emph{SIAM Rev.}, 43(1): pp.  89--112, 2001.

\bibitem{Gough_GSL}
B.~Gough.
\newblock GNU Scientific Library Reference Manual - Third Edition.
\newblock Network Theory Ltd., 3rd edition, 2009.
\newblock ISBN 0954612078, 9780954612078.

\bibitem{hesthaven2007nodal}
J.~S. Hesthaven and T.~Warburton.
\newblock Nodal discontinuous Galerkin methods: algorithms, analysis, and
  applications.
\newblock Springer Science \& Business Media, 2007.

\bibitem{jaynes1957info}
E.~T. Jaynes.
\newblock Information theory and statistical mechanics.
\newblock \emph{Phys. Rev.}, 106(4):page 620, 1957.

\bibitem{Kanevsky2007}
A.~Kanevsky, M.~H. Carpenter, D.~Gottlieb, and J.~S. Hesthaven.
\newblock {Application of implicit–explicit high order Runge–Kutta methods
  to discontinuous-Galerkin schemes}.
\newblock \emph{Journal of Computational Physics}, 225(2): pp.  1753--1781,
  2007.
\newblock ISSN 00219991.

\bibitem{Ker76}
D.~Kershaw.
\newblock Flux limiting nature's own way: A new method for numerical solution
  of the transport equation.
\newblock Technical report, LLNL Report UCRL-78378, 1976.

\bibitem{Ketcheson2008}
D.~I. Ketcheson.
\newblock {Highly Efficient Strong Stability-Preserving Runge–Kutta Methods
  with Low-Storage Implementations}.
\newblock \emph{SIAM Journal on Scientific Computing}, 30(4): pp.  2113--2136,
  2008.
\newblock ISSN 1064-8275.

\bibitem{Ketcheson}
D.~I. Ketcheson, S.~Gottlieb, and C.~B. Macdonald.
\newblock {Strong stability preserving two-step Runge-Kutta methods}.
\newblock \emph{SIAM Journal on Numerical {\ldots}},  pp.  1--24, 2011.

\bibitem{Krivodonova2007}
L.~Krivodonova.
\newblock Limiters for high-order discontinuous galerkin methods.
\newblock \emph{Journal of Computational Physics}, 226(1): pp.  879 -- 896,
  2007.
\newblock ISSN 0021-9991.

\bibitem{Krivodonova2004}
L.~Krivodonova, J.~Xin, J.-F. Remacle, N.~Chevaugeon, and J.~Flaherty.
\newblock {Shock detection and limiting with discontinuous Galerkin methods for
  hyperbolic conservation laws}.
\newblock \emph{Appl. Numer. Math.}, 48(3-4): pp.  323--338, 2004.

\bibitem{Larsen1991}
E.~W. Larsen and G.~C. Pomraning.
\newblock {The PN Theory as an Asymptotic Limit of Transport Theory in Planar
  Geometry —I: Analysis}.
\newblock \emph{Nuclear Science and Engineering}, 109(1): pp.  49--75, 1991.

\bibitem{LeVeque2002}
R.~LeVeque.
\newblock {Finite volume methods for hyperbolic problems}.
\newblock 2002.

\bibitem{Lev84}
C.~D. Levermore.
\newblock Relating {E}ddington factors to flux limiters.
\newblock \emph{J. Quant. Spectrosc. Radiat. Transfer}, 31: pp.  149--160,
  1984.

\bibitem{levermore1996moment}
---.
\newblock Moment closure hierarchies for kinetic theories.
\newblock \emph{J. Stat. Phys.}, 83(5--6): pp.  1021--1065, 1996.

\bibitem{levermore2009boundary}
---.
\newblock {Boundary conditions for moment closures}.
\newblock \emph{Institute for Pure and Applied Mathematics University of
  California, Los Angeles, CA on May}, 27, 2009.

\bibitem{li2011central}
F.~Li, L.~Xu, and S.~Yakovlev.
\newblock {C}entral discontinuous {G}alerkin methods for ideal {MHD} equations
  with the exactly divergence-free magnetic field.
\newblock \emph{J. Comput. Phys.}, 230(12): pp.  4828--4847, 2011.

\bibitem{liangdiscontinuous}
C.~Liang, F.~Ham, and E.~Johnsen.
\newblock {D}iscontinuous {G}alerkin method with {WENO} limiter for flows with
  discontinuity.
\newblock Technical report, Center for Turbulence Research, 2009.

\bibitem{liu_xu_gas_kinetic}
H.~Liu, , and K.~Xu.
\newblock A gas-kinetic discontinuous {G}alerkin method for viscous flow
  equations.
\newblock \emph{J. Mech. Sci. Technol.}, 21(9): pp.  1344--1351, 2007.

\bibitem{liu2007runge}
H.~Liu and K.~Xu.
\newblock {A} {R}unge-{K}utta discontinuous {G}alerkin method for viscous flow
  equations.
\newblock \emph{J. Comput. Phys.}, 224(2): pp.  1223--1242, 2007.

\bibitem{Monreal}
P.~Monreal.
\newblock Moment Realizability and Kershaw Closures in Radiative Transfer.
\newblock Ph.D. thesis, TU Aachen, 2012.

\bibitem{MullerRuggeri}
I.~M{\"u}ller and T.~Ruggeri.
\newblock Rational Extended Thermodynamics.
\newblock Springer-Verlag, New York, 2nd edition, 1993.

\bibitem{Olbrant12}
E.~Olbrant, C.~D. Hauck, and M.~Frank.
\newblock A realizability-preserving discontinuous galerkin method for the
  ${M}_1$ model of radiative transfer.
\newblock \emph{J. Comput. Phys.}, 231(17): pp.  5612 -- 5639, 2012.

\bibitem{Ore55}
A.~Ore.
\newblock Entropy of radiation.
\newblock \emph{Phys. Rev.}, 98:page 887, 1955.

\bibitem{pomraning1964variational}
G.~C. Pomraning.
\newblock {Variational boundary conditions for the spherical harmonics
  approximation to the neutron transport equation}.
\newblock \emph{Annals of Physics}, 27(2): pp.  193--215, 1964.

\bibitem{Qiu2005a}
J.~Qiu and C.~Shu.
\newblock {A comparison of troubled cell indicators for Runge-Kutta
  discontinuous Galerkin methods using WENO Limiters}.
\newblock \emph{SIAM J. Sci. Comput.}, 27: pp.  995--1013, 2005.

\bibitem{Radice2013648}
D.~Radice, E.~Abdikamalov, L.~Rezzolla, and C.~D. Ott.
\newblock {A new spherical harmonics scheme for multi-dimensional radiation
  transport I. Static matter configurations}.
\newblock \emph{J. Comput. Phys.}, 242(0): pp.  648--669, 2013.

\bibitem{Rampp2002}
M.~Rampp and H.~T. Janka.
\newblock {Radiation hydrodynamics with neutrinos: Variable Eddington factor
  method for core-collapse supernova simulations}.
\newblock \emph{Astron. Astrophys.}, 396:page~42, 2002.

\bibitem{reedhill}
W.~Reed and T.~Hill.
\newblock Triangular mesh methods for neutron transport equation.
\newblock \emph{Los Alamos Scientific Laboratory Report}, LA-UR-73-479, 1973.

\bibitem{Ros54}
P.~Rosen.
\newblock Entropy of radiation.
\newblock \emph{Phys. Rev.}, 96:page 555, 1954.

\bibitem{Rulko1991}
R.~P. Rulko, E.~W. Larsen, and G.~C. Pomraning.
\newblock {The PN Theory as an Asymptotic Limit of Transport Theory in Planar
  Geometry —II: Numerical Results}.
\newblock \emph{Nuclear Science and Engineering}, 109(1): pp.  76--85, 1991.

\bibitem{Schneider2016}
F.~Schneider.
\newblock {Moment models in radiation transport equations}.
\newblock Dr. Hut Verlag, 2016.

\bibitem{Science2013}
F.~Science, Y.~Kanno, T.~Harada, and T.~Hanawa.
\newblock {Kinetic Scheme for Solving the M1 Model of Radiative Transfer}.
\newblock \emph{Publ. Astron. Soc. Jpn.}, 65(2007):page~72, 2013.

\bibitem{CWShu}
C.~Shu.
\newblock {TVD} time discretizations.
\newblock \emph{SIAM J. Sci. Stat. Comput.}, 9: pp.  1073--1084, 1988.

\bibitem{Shu1998}
C.-W. Shu.
\newblock {Essentially non-oscillatory and weighted essentially non-oscillatory
  schemes for hyperbolic conservation laws}.
\newblock 1998.

\bibitem{Shu1988}
C.-W. Shu and S.~Osher.
\newblock Efficient implementation of essentially non-oscillatory
  shock-capturing schemes.
\newblock \emph{J. Comput. Phys.}, 77(2): pp.  439--471, 1988.

\bibitem{smit1997hyperbolicity}
J.~M. Smit, J.~Cernohorsky, and C.~P. Dullemond.
\newblock {Hyperbolicity and critical points in two-moment approximate
  radiative transfer.}
\newblock \emph{Astron. Astrophys.}, 325: pp.  203--211, 1997.

\bibitem{Struchtrup2000}
H.~Struchtrup.
\newblock {Kinetic schemes and boundary conditions for moment equations}.
\newblock \emph{Zeitschrift f{\"{u}}r angewandte Mathematik und Physik},
  51(3):page 346, 2000.
\newblock ISSN 00442275.
\newblock \doi{10.1007/s000330050002}.

\bibitem{Sumiyoshi2012}
K.~Sumiyoshi and S.~Yamada.
\newblock {Neutrino Transfer in Three Dimensions for Core-Collapse Supernovae.
  I. Static Configurations}.
\newblock \emph{Astrophys. J. Suppl. Ser.}, 199(1):page~17, 2012.

\bibitem{toro2009riemann}
E.~Toro.
\newblock Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
  Introduction.
\newblock Springer, 2009.
\newblock ISBN 9783540252023.

\bibitem{Wang2012}
C.~Wang, X.~Zhang, C.~W. Shu, and J.~Ning.
\newblock {Robust high order discontinuous Galerkin schemes for two-dimensional
  gaseous detonations}.
\newblock \emph{J. Comput. Phys.}, 231(2): pp.  653--665, 2012.

\bibitem{Xing2013}
Y.~Xing and X.~Zhang.
\newblock {Positivity-preserving well-balanced discontinuous Galerkin methods
  for the shallow water equations on unstructured triangular meshes}.
\newblock \emph{J. Sci. Comput.},  pp.  1--35, 2013.

\bibitem{Zhang2010}
X.~Zhang and C.~W. Shu.
\newblock {On positivity-preserving high order discontinuous Galerkin schemes
  for compressible Euler equations on rectangular meshes}.
\newblock \emph{J. Comput. Phys.}, 229(23): pp.  8918--8934, 2010.

\bibitem{Zhang2012a}
X.~Zhang and C.-W. Shu.
\newblock {A minimum entropy principle of high order schemes for gas dynamics
  equations}.
\newblock \emph{Numer. Math.}, 121(3): pp.  545--563, 2012.

\bibitem{ZhangShu12}
X.~Zhang, Y.~Xia, and C.-W. Shu.
\newblock Maximum-principle-satisfying and positivity-preserving high order
  discontinuous {G}alerkin schemes for conservation laws on triangular meshes.
\newblock \emph{J. Sci. Comput.}, 50(1): pp.  29--62, 2012.

\end{thebibliography}
