subroutine monomial_value_2d ( n, ex, ey, x, y, v ) !*****************************************************************************80 ! !! monomial_value_2d() evaluates a monomial in x and y. ! ! Discussion: ! ! This routine evaluates a monomial of the form ! ! x^ex * y^ey ! ! The combination 0.0^0 is encountered is treated as 1.0. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, integer EX, EY, the exponents. ! ! Input, real ( kind = rk ) X(N), Y(N), the point coordinates. ! ! Output, real ( kind = rk ) V(N), the monomial values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer ex integer ey real ( kind = rk ) v(n) real ( kind = rk ) x(n) real ( kind = rk ) y(n) v(1:n) = 1.0D+00 if ( 0 /= ex ) then v(1:n) = v(1:n) * x(1:n) ** ex end if if ( 0 /= ey ) then v(1:n) = v(1:n) * y(1:n) ** ey end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine triangle_unit_monomial ( ex, ey, value ) !*****************************************************************************80 ! !! triangle_unit_monomial integrates a monomial over the unit triangle. ! ! Discussion: ! ! This routine integrates a monomial of the form ! ! x^ex y^ey ! ! where the exponents are nonnegative integers. Note that ! if the combination 0^0 is encountered, it should be treated ! as 1. ! ! Integral ( over unit triangle ) x^m y^n dx dy = m! * n! / ( m + n + 2 )! ! ! The integration region is: ! ! 0 <= X ! 0 <= Y ! X + Y <= 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer EX, EY, the exponents. ! ! Output, real ( kind = rk ) VALUE, the integral of the monomial. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ex integer ey integer i integer k real ( kind = rk ) value value = 1.0D+00 k = ex do i = 1, ey k = k + 1 value = value * real ( i, kind = rk ) / real ( k, kind = rk ) end do k = k + 1 value = value / real ( k, kind = rk ) k = k + 1 value = value / real ( k, kind = rk ) return end function triangle_unit_volume ( ) !*****************************************************************************80 ! !! TRIANGLE_UNIT_VOLUME: volume of a unit triangle. ! ! Discussion: ! ! The integration region is: ! ! 0 <= X ! 0 <= Y ! X + Y <= 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 March 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) TRIANGLE_UNIT_VOLUME, the volume. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) triangle_unit_volume triangle_unit_volume = 0.5D+00 return end subroutine twb_rule_n ( strength, n ) !*****************************************************************************80 ! !! twb_rule_n returns the order of a TWB rule of given strength. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Mark Taylor, Beth Wingate, Len Bos, ! Several new quadrature formulas for polynomial integration in the triangle, ! http://arxiv.org/abs/math/0501496v2, ! 08 February 2007. ! ! Parameters: ! ! Input, integer STRENGTH, the desired strength. ! 1 <= STRENGTH. ! ! Output, integer N, the order of the rule, if it exists. ! Otherwise, the value -1 is returned. ! implicit none integer n integer strength if ( strength == 1 ) then n = 1 else if ( strength == 2 ) then n = 3 else if ( strength == 4 ) then n = 6 else if ( strength == 5 ) then n = 10 else if ( strength == 7 ) then n = 15 else if ( strength == 9 ) then n = 21 else if ( strength == 11 ) then n = 28 else if ( strength == 13 ) then n = 36 else if ( strength == 14 ) then n = 45 else if ( strength == 16 ) then n = 55 else if ( strength == 18 ) then n = 66 else if ( strength == 20 ) then n = 78 else if ( strength == 21 ) then n = 91 else if ( strength == 23 ) then n = 105 else if ( strength == 25 ) then n = 120 else n = -1 end if return end subroutine twb_rule_w ( strength, w ) !*****************************************************************************80 ! !! twb_rule_w returns the weight of a TWB rule of given strength. ! ! Discussion: ! ! The weights are normalized so that they sum to 1/2, the area of the ! triangle (0,0) -- (1,0) -- (1,0). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Mark Taylor, Beth Wingate, Len Bos, ! Several new quadrature formulas for polynomial integration in the triangle, ! http://arxiv.org/abs/math/0501496v2, ! 08 February 2007. ! ! Parameters: ! ! Input, integer STRENGTH, the desired strength. ! 1 <= STRENGTH. ! ! Output, real ( kind = rk ) W(N), the weight of the rule, if it exists. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer strength real ( kind = rk ) w(*) call twb_rule_n ( strength, n ) if ( n < 1 ) then return end if if ( strength == 1 ) then w(1:n) = (/ & 2.0D+00 /) else if ( strength == 2 ) then w(1:n) = (/ & 0.6666666666667D+00, & 0.6666666666667D+00, & 0.6666666666667D+00 /) else if ( strength == 4 ) then w(1:n) = (/ & 0.2199034873106D+00, & 0.2199034873106D+00, & 0.2199034873106D+00, & 0.4467631793560D+00, & 0.4467631793560D+00, & 0.4467631793560D+00 /) else if ( strength == 5 ) then w(1:n) = (/ & 0.0262712099504D+00, & 0.0262716612068D+00, & 0.0274163947600D+00, & 0.2348383865823D+00, & 0.2348412238268D+00, & 0.2480251793114D+00, & 0.2480304922521D+00, & 0.2518604605529D+00, & 0.2518660533658D+00, & 0.4505789381914D+00 /) else if ( strength == 7 ) then w(1:n) = (/ & 0.0102558174092D+00, & 0.0102558174092D+00, & 0.0102558174092D+00, & 0.1116047046647D+00, & 0.1116047046647D+00, & 0.1116047046647D+00, & 0.1116047046647D+00, & 0.1116047046647D+00, & 0.1116047046647D+00, & 0.1679775595335D+00, & 0.1679775595335D+00, & 0.1679775595335D+00, & 0.2652238803946D+00, & 0.2652238803946D+00, & 0.2652238803946D+00 /) else if ( strength == 9 ) then w(1:n) = (/ & 0.0519871420646D+00, & 0.0519871420646D+00, & 0.0519871420646D+00, & 0.0707034101784D+00, & 0.0707034101784D+00, & 0.0707034101784D+00, & 0.0707034101784D+00, & 0.0707034101784D+00, & 0.0707034101784D+00, & 0.0909390760952D+00, & 0.0909390760952D+00, & 0.0909390760952D+00, & 0.0909390760952D+00, & 0.0909390760952D+00, & 0.0909390760952D+00, & 0.1032344051380D+00, & 0.1032344051380D+00, & 0.1032344051380D+00, & 0.1881601469167D+00, & 0.1881601469167D+00, & 0.1881601469167D+00 /) else if ( strength == 11 ) then w(1:n) = (/ & 0.0114082494033D+00, & 0.0114082494033D+00, & 0.0132691285720D+00, & 0.0132691285720D+00, & 0.0155865773350D+00, & 0.0408274780428D+00, & 0.0408274780429D+00, & 0.0579849665116D+00, & 0.0579849665116D+00, & 0.0601385247663D+00, & 0.0601385247663D+00, & 0.0625273888433D+00, & 0.0625273888433D+00, & 0.0639684321504D+00, & 0.0639684321504D+00, & 0.0661325872161D+00, & 0.0668503236820D+00, & 0.0668503236821D+00, & 0.0686904305977D+00, & 0.0686904305977D+00, & 0.1002717543859D+00, & 0.1143136784099D+00, & 0.1143136784099D+00, & 0.1223648146752D+00, & 0.1223648146752D+00, & 0.1394422334178D+00, & 0.1394422334178D+00, & 0.1744377829182D+00 /) else if ( strength == 13 ) then w(1:n) = (/ & 0.0166240998757D+00, & 0.0166811699778D+00, & 0.0166830569067D+00, & 0.0175680870083D+00, & 0.0184474661845D+00, & 0.0197942410188D+00, & 0.0203540395855D+00, & 0.0206852863940D+00, & 0.0208271366086D+00, & 0.0317819778279D+00, & 0.0320472035241D+00, & 0.0320607681146D+00, & 0.0430765959183D+00, & 0.0438473415339D+00, & 0.0439209672733D+00, & 0.0479951923691D+00, & 0.0483806260733D+00, & 0.0484867423375D+00, & 0.0556964488024D+00, & 0.0561026364356D+00, & 0.0565190123693D+00, & 0.0689289890670D+00, & 0.0717213336089D+00, & 0.0727453920976D+00, & 0.0788807336737D+00, & 0.0810114345512D+00, & 0.0825725299055D+00, & 0.0842044567330D+00, & 0.0843585533305D+00, & 0.0851969868488D+00, & 0.0902845328052D+00, & 0.0914283143485D+00, & 0.0916279065409D+00, & 0.1025573374896D+00, & 0.1033159661413D+00, & 0.1035854367193D+00 /) else if ( strength == 14 ) then w(1:n) = (/ & 0.0010616711990D+00, & 0.0010616711990D+00, & 0.0010616711990D+00, & 0.0131460236101D+00, & 0.0131460236101D+00, & 0.0131460236101D+00, & 0.0131460236101D+00, & 0.0131460236101D+00, & 0.0131460236101D+00, & 0.0242881926949D+00, & 0.0242881926949D+00, & 0.0242881926949D+00, & 0.0242881926949D+00, & 0.0242881926949D+00, & 0.0242881926949D+00, & 0.0316799866332D+00, & 0.0316799866332D+00, & 0.0316799866332D+00, & 0.0316799866332D+00, & 0.0316799866332D+00, & 0.0316799866332D+00, & 0.0349317947036D+00, & 0.0349317947036D+00, & 0.0349317947036D+00, & 0.0383664533945D+00, & 0.0383664533945D+00, & 0.0383664533945D+00, & 0.0578369491210D+00, & 0.0578369491210D+00, & 0.0578369491210D+00, & 0.0578369491210D+00, & 0.0578369491210D+00, & 0.0578369491210D+00, & 0.0725821687394D+00, & 0.0725821687394D+00, & 0.0725821687394D+00, & 0.0725821687394D+00, & 0.0725821687394D+00, & 0.0725821687394D+00, & 0.0897856524107D+00, & 0.0897856524107D+00, & 0.0897856524107D+00, & 0.1034544533617D+00, & 0.1034544533617D+00, & 0.1034544533617D+00 /) else if ( strength == 16 ) then w(1:n) = (/ & 0.0006202599851D+00, & 0.0006315174712D+00, & 0.0007086601559D+00, & 0.0055163716168D+00, & 0.0062692407656D+00, & 0.0078531408826D+00, & 0.0094551483864D+00, & 0.0097824511271D+00, & 0.0099861643489D+00, & 0.0137553818816D+00, & 0.0140979178040D+00, & 0.0149646864337D+00, & 0.0156097503612D+00, & 0.0157683693348D+00, & 0.0175794546383D+00, & 0.0204113840270D+00, & 0.0209562878616D+00, & 0.0210713412998D+00, & 0.0217646760202D+00, & 0.0222288408699D+00, & 0.0224186693682D+00, & 0.0230122616993D+00, & 0.0236813902500D+00, & 0.0257464643368D+00, & 0.0257956801608D+00, & 0.0258072327610D+00, & 0.0260343232059D+00, & 0.0265768141609D+00, & 0.0265784761831D+00, & 0.0267532329238D+00, & 0.0375787806641D+00, & 0.0383065894195D+00, & 0.0384849695025D+00, & 0.0389619825852D+00, & 0.0394604111547D+00, & 0.0412364778098D+00, & 0.0512872438483D+00, & 0.0516405641935D+00, & 0.0518230042269D+00, & 0.0528527988181D+00, & 0.0538505573027D+00, & 0.0541895329319D+00, & 0.0584737146444D+00, & 0.0592863168363D+00, & 0.0594358276749D+00, & 0.0631800255863D+00, & 0.0632926845153D+00, & 0.0640707361772D+00, & 0.0812040595918D+00, & 0.0814437513530D+00, & 0.0814679201241D+00, & 0.0815050548084D+00, & 0.0815164664939D+00, & 0.0816931059623D+00, & 0.0923218334531D+00 /) else if ( strength == 18 ) then w(1:n) = (/ & 0.0025165756986D+00, & 0.0025273452007D+00, & 0.0033269295333D+00, & 0.0081503492125D+00, & 0.0086135525742D+00, & 0.0087786746179D+00, & 0.0097099585562D+00, & 0.0102466211915D+00, & 0.0108397688341D+00, & 0.0129385390176D+00, & 0.0136339823583D+00, & 0.0138477328147D+00, & 0.0139421540105D+00, & 0.0144121399968D+00, & 0.0153703455534D+00, & 0.0162489802253D+00, & 0.0169718304280D+00, & 0.0170088532421D+00, & 0.0170953520675D+00, & 0.0173888854559D+00, & 0.0174543962439D+00, & 0.0178406757287D+00, & 0.0178446863879D+00, & 0.0179046337552D+00, & 0.0181259756201D+00, & 0.0184784838882D+00, & 0.0185793564371D+00, & 0.0203217151777D+00, & 0.0213771661809D+00, & 0.0231916854098D+00, & 0.0274426710859D+00, & 0.0290301922340D+00, & 0.0294522738505D+00, & 0.0299436251629D+00, & 0.0307026948119D+00, & 0.0325263365863D+00, & 0.0327884208506D+00, & 0.0331234675192D+00, & 0.0346167526875D+00, & 0.0347081373976D+00, & 0.0347372049404D+00, & 0.0348528762454D+00, & 0.0348601561186D+00, & 0.0355471569975D+00, & 0.0360182996383D+00, & 0.0362926285843D+00, & 0.0381897702083D+00, & 0.0392252800118D+00, & 0.0482710125888D+00, & 0.0489912121566D+00, & 0.0497220833872D+00, & 0.0507065736986D+00, & 0.0509771994043D+00, & 0.0521360063667D+00, & 0.0523460874925D+00, & 0.0524440683552D+00, & 0.0527459644823D+00, & 0.0529449063728D+00, & 0.0542395594501D+00, & 0.0543470203419D+00, & 0.0547100548639D+00, & 0.0557288345913D+00, & 0.0577734264233D+00, & 0.0585393781623D+00, & 0.0609039250680D+00, & 0.0637273964449D+00 /) else if ( strength == 20 ) then w(1:n) = (/ & 0.0021744545399D+00, & 0.0028987135265D+00, & 0.0030846029337D+00, & 0.0034401633104D+00, & 0.0041898472012D+00, & 0.0044738051498D+00, & 0.0047054420814D+00, & 0.0048867935750D+00, & 0.0051927643369D+00, & 0.0074073058981D+00, & 0.0079755410301D+00, & 0.0083550522910D+00, & 0.0096166660864D+00, & 0.0096318257850D+00, & 0.0098577460758D+00, & 0.0102657880301D+00, & 0.0103188103111D+00, & 0.0106291001630D+00, & 0.0106881306895D+00, & 0.0106969021010D+00, & 0.0109026461714D+00, & 0.0109899783575D+00, & 0.0113423055229D+00, & 0.0120535642930D+00, & 0.0139619193821D+00, & 0.0141147991536D+00, & 0.0141930347046D+00, & 0.0144212676268D+00, & 0.0144704346855D+00, & 0.0144949769872D+00, & 0.0145386775694D+00, & 0.0145964190926D+00, & 0.0147314578466D+00, & 0.0167463963304D+00, & 0.0168955500458D+00, & 0.0169422662884D+00, & 0.0173070172095D+00, & 0.0174524546493D+00, & 0.0177217222159D+00, & 0.0282824024023D+00, & 0.0284996712488D+00, & 0.0285005646539D+00, & 0.0300647223478D+00, & 0.0302031277082D+00, & 0.0303987136077D+00, & 0.0305668796074D+00, & 0.0306067413002D+00, & 0.0309330068201D+00, & 0.0309773820835D+00, & 0.0313146250545D+00, & 0.0313573493392D+00, & 0.0314320469287D+00, & 0.0315182143894D+00, & 0.0324248137985D+00, & 0.0347512152386D+00, & 0.0350393454927D+00, & 0.0350717420310D+00, & 0.0352129215334D+00, & 0.0352615504981D+00, & 0.0366403220343D+00, & 0.0367733107670D+00, & 0.0371675662937D+00, & 0.0373371571606D+00, & 0.0403973346588D+00, & 0.0413580040638D+00, & 0.0421957791870D+00, & 0.0495451004037D+00, & 0.0500419261141D+00, & 0.0505794587115D+00, & 0.0520037210188D+00, & 0.0521533567886D+00, & 0.0524899152358D+00, & 0.0599159762516D+00, & 0.0599609997426D+00, & 0.0599915272129D+00, & 0.0634133183449D+00, & 0.0635311861108D+00, & 0.0637206605672D+00 /) else if ( strength == 21 ) then w(1:n) = (/ & 0.0006704436439D+00, & 0.0006704436439D+00, & 0.0006704436439D+00, & 0.0045472608074D+00, & 0.0045472608074D+00, & 0.0045472608074D+00, & 0.0045472608074D+00, & 0.0045472608074D+00, & 0.0045472608074D+00, & 0.0052077585320D+00, & 0.0052077585320D+00, & 0.0052077585320D+00, & 0.0052077585320D+00, & 0.0052077585320D+00, & 0.0052077585320D+00, & 0.0065435432887D+00, & 0.0065435432887D+00, & 0.0065435432887D+00, & 0.0092737841533D+00, & 0.0092737841533D+00, & 0.0092737841533D+00, & 0.0092737841533D+00, & 0.0092737841533D+00, & 0.0092737841533D+00, & 0.0095937782623D+00, & 0.0095937782623D+00, & 0.0095937782623D+00, & 0.0095937782623D+00, & 0.0095937782623D+00, & 0.0095937782623D+00, & 0.0114247809167D+00, & 0.0114247809167D+00, & 0.0114247809167D+00, & 0.0114247809167D+00, & 0.0114247809167D+00, & 0.0114247809167D+00, & 0.0117216964174D+00, & 0.0117216964174D+00, & 0.0117216964174D+00, & 0.0117216964174D+00, & 0.0117216964174D+00, & 0.0117216964174D+00, & 0.0188197155232D+00, & 0.0188197155232D+00, & 0.0188197155232D+00, & 0.0188197155232D+00, & 0.0188197155232D+00, & 0.0188197155232D+00, & 0.0235260980271D+00, & 0.0235260980271D+00, & 0.0235260980271D+00, & 0.0235571466151D+00, & 0.0235571466151D+00, & 0.0235571466151D+00, & 0.0268246207430D+00, & 0.0268246207430D+00, & 0.0268246207430D+00, & 0.0268246207430D+00, & 0.0268246207430D+00, & 0.0268246207430D+00, & 0.0314289776779D+00, & 0.0314289776779D+00, & 0.0314289776779D+00, & 0.0314289776779D+00, & 0.0314289776779D+00, & 0.0314289776779D+00, & 0.0337196192159D+00, & 0.0337196192159D+00, & 0.0337196192159D+00, & 0.0337196192159D+00, & 0.0337196192159D+00, & 0.0337196192159D+00, & 0.0427745294213D+00, & 0.0427745294213D+00, & 0.0427745294213D+00, & 0.0427745294213D+00, & 0.0427745294213D+00, & 0.0427745294213D+00, & 0.0441138932737D+00, & 0.0441138932737D+00, & 0.0441138932737D+00, & 0.0461469594684D+00, & 0.0461469594684D+00, & 0.0461469594684D+00, & 0.0461469594684D+00, & 0.0461469594684D+00, & 0.0461469594684D+00, & 0.0469152468624D+00, & 0.0469152468624D+00, & 0.0469152468624D+00, & 0.0551199980347D+00 /) else if ( strength == 23 ) then w(1:n) = (/ & 0.0006438298261D+00, & 0.0006438413076D+00, & 0.0010134735710D+00, & 0.0010134752576D+00, & 0.0019679929935D+00, & 0.0033467313784D+00, & 0.0033467339208D+00, & 0.0042873323375D+00, & 0.0042873459885D+00, & 0.0043003801372D+00, & 0.0043003849098D+00, & 0.0056934629205D+00, & 0.0056934640134D+00, & 0.0061643868015D+00, & 0.0061644756418D+00, & 0.0062014513591D+00, & 0.0062014531952D+00, & 0.0069636330294D+00, & 0.0069636331842D+00, & 0.0075066257720D+00, & 0.0075066264565D+00, & 0.0079074768339D+00, & 0.0079074772485D+00, & 0.0080353344623D+00, & 0.0087963441074D+00, & 0.0087963448112D+00, & 0.0091304195716D+00, & 0.0091304213611D+00, & 0.0092821748751D+00, & 0.0092821815662D+00, & 0.0094499806178D+00, & 0.0094627468484D+00, & 0.0094627485294D+00, & 0.0095555772285D+00, & 0.0095555792843D+00, & 0.0096138842488D+00, & 0.0096138846826D+00, & 0.0099991524212D+00, & 0.0099991551850D+00, & 0.0100301319277D+00, & 0.0100301346636D+00, & 0.0124936676185D+00, & 0.0124936726125D+00, & 0.0140197309137D+00, & 0.0143336216896D+00, & 0.0143336272125D+00, & 0.0153604142740D+00, & 0.0153604183425D+00, & 0.0184523825614D+00, & 0.0184523863146D+00, & 0.0195833983573D+00, & 0.0195834019994D+00, & 0.0197632751342D+00, & 0.0197632766677D+00, & 0.0198806391019D+00, & 0.0198806485776D+00, & 0.0207181838484D+00, & 0.0207181934893D+00, & 0.0208943071440D+00, & 0.0208943251956D+00, & 0.0214864573885D+00, & 0.0214864586007D+00, & 0.0222218133036D+00, & 0.0222218160203D+00, & 0.0223345305455D+00, & 0.0223345378739D+00, & 0.0224758924946D+00, & 0.0224758980440D+00, & 0.0229701395845D+00, & 0.0229703394438D+00, & 0.0232798376102D+00, & 0.0232798427506D+00, & 0.0269483199647D+00, & 0.0269483307107D+00, & 0.0280438758010D+00, & 0.0280438764607D+00, & 0.0287526270172D+00, & 0.0287526387271D+00, & 0.0298980829063D+00, & 0.0298980922759D+00, & 0.0309004358516D+00, & 0.0309004385956D+00, & 0.0314031017088D+00, & 0.0314031073955D+00, & 0.0319191553024D+00, & 0.0319191668378D+00, & 0.0321429924062D+00, & 0.0330395601388D+00, & 0.0330395631829D+00, & 0.0356169095589D+00, & 0.0356169276054D+00, & 0.0365741189998D+00, & 0.0365741515204D+00, & 0.0365977646990D+00, & 0.0365978053889D+00, & 0.0369945680114D+00, & 0.0369945775059D+00, & 0.0374053623787D+00, & 0.0375550258317D+00, & 0.0375550312530D+00, & 0.0388887693486D+00, & 0.0388887708342D+00, & 0.0392705643548D+00, & 0.0392705802517D+00, & 0.0398766879831D+00 /) else if ( strength == 25 ) then w(1:n) = (/ & 0.0014873417859D+00, & 0.0014889035262D+00, & 0.0015005944380D+00, & 0.0015059208313D+00, & 0.0015318868715D+00, & 0.0023032634487D+00, & 0.0023649067042D+00, & 0.0028751143611D+00, & 0.0029862488735D+00, & 0.0030384162737D+00, & 0.0032092459688D+00, & 0.0037029598435D+00, & 0.0037407186035D+00, & 0.0038452543223D+00, & 0.0038670778668D+00, & 0.0039192555178D+00, & 0.0039573282688D+00, & 0.0044032251724D+00, & 0.0045907108173D+00, & 0.0047023669435D+00, & 0.0050014843818D+00, & 0.0052387830156D+00, & 0.0054422104092D+00, & 0.0056931248912D+00, & 0.0059107422989D+00, & 0.0059687967687D+00, & 0.0067262190287D+00, & 0.0068307848624D+00, & 0.0069531259112D+00, & 0.0072460270642D+00, & 0.0072728189613D+00, & 0.0073008930847D+00, & 0.0073604666776D+00, & 0.0074119923255D+00, & 0.0074892214336D+00, & 0.0078604067260D+00, & 0.0078621726423D+00, & 0.0080506361066D+00, & 0.0081442860473D+00, & 0.0081478804152D+00, & 0.0092444146612D+00, & 0.0094674635165D+00, & 0.0097132210137D+00, & 0.0099753581151D+00, & 0.0103367803673D+00, & 0.0112263277166D+00, & 0.0114309118745D+00, & 0.0115550567487D+00, & 0.0135575856957D+00, & 0.0135984962900D+00, & 0.0137754813837D+00, & 0.0137961015942D+00, & 0.0138408839904D+00, & 0.0140634019977D+00, & 0.0140991451009D+00, & 0.0142004111991D+00, & 0.0144518424517D+00, & 0.0150245979639D+00, & 0.0152817804122D+00, & 0.0155550724169D+00, & 0.0164570886000D+00, & 0.0165275759573D+00, & 0.0166847554451D+00, & 0.0167409312985D+00, & 0.0168674663361D+00, & 0.0168882230165D+00, & 0.0172087112691D+00, & 0.0174681068264D+00, & 0.0176663899614D+00, & 0.0182967621475D+00, & 0.0183576852459D+00, & 0.0186392569521D+00, & 0.0189781060590D+00, & 0.0191847922578D+00, & 0.0194080442044D+00, & 0.0194720072193D+00, & 0.0200855080495D+00, & 0.0201673909332D+00, & 0.0221742162761D+00, & 0.0229702440508D+00, & 0.0233465117399D+00, & 0.0234883135338D+00, & 0.0240682099018D+00, & 0.0240910792953D+00, & 0.0245677049481D+00, & 0.0246536315719D+00, & 0.0246756530052D+00, & 0.0249704602710D+00, & 0.0250026544082D+00, & 0.0250490869426D+00, & 0.0250936250125D+00, & 0.0251482076226D+00, & 0.0255010290447D+00, & 0.0256544511979D+00, & 0.0257974750630D+00, & 0.0270007753993D+00, & 0.0274431536844D+00, & 0.0277072401488D+00, & 0.0278284415364D+00, & 0.0287207381105D+00, & 0.0288826834956D+00, & 0.0293302729759D+00, & 0.0318902879557D+00, & 0.0319083660286D+00, & 0.0320938960329D+00, & 0.0321618608780D+00, & 0.0322424127534D+00, & 0.0327072446421D+00, & 0.0329946316695D+00, & 0.0331828096025D+00, & 0.0334857162651D+00, & 0.0335468472792D+00, & 0.0337049042988D+00, & 0.0340361462767D+00, & 0.0342465235323D+00, & 0.0345528817251D+00, & 0.0356782875703D+00, & 0.0364656225016D+00, & 0.0365172708706D+00, & 0.0371924811018D+00 /) end if w(1:n) = w(1:n) / 4.0D+00 return end subroutine twb_rule_x ( strength, x ) !*****************************************************************************80 ! !! twb_rule_x returns the x abscissas of a TWB rule of given strength. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Mark Taylor, Beth Wingate, Len Bos, ! Several new quadrature formulas for polynomial integration in the triangle, ! http://arxiv.org/abs/math/0501496v2, ! 08 February 2007. ! ! Parameters: ! ! Input, integer STRENGTH, the desired strength. ! 1 <= STRENGTH. ! ! Output, real ( kind = rk ) X(N), the x abscissas of the rule, if it exists. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer strength real ( kind = rk ) x(*) call twb_rule_n ( strength, n ) if ( n < 1 ) then return end if if ( strength == 1 ) then x(1:n) = (/ & 0.3333333333333D+00 /) elseif ( strength == 2 ) then x(1:n) = (/ & 0.1666666666667D+00, & 0.6666666666667D+00, & 0.1666666666667D+00 /) elseif ( strength == 4 ) then x(1:n) = (/ & 0.0915762135098D+00, & 0.8168475729805D+00, & 0.0915762135098D+00, & 0.1081030181681D+00, & 0.4459484909160D+00, & 0.4459484909160D+00 /) elseif ( strength == 5 ) then x(1:n) = (/ & 0.0000000000000D+00, & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.2673273531185D+00, & 0.6728175529461D+00, & 0.0649236350054D+00, & 0.6716498539042D+00, & 0.0654032456800D+00, & 0.2693767069140D+00, & 0.3386738503896D+00 /) elseif ( strength == 7 ) then x(1:n) = (/ & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0000000000000D+00, & 0.7839656651012D+00, & 0.1738960507345D+00, & 0.1738960507345D+00, & 0.0421382841642D+00, & 0.7839656651012D+00, & 0.0421382841642D+00, & 0.4743880861752D+00, & 0.4743880861752D+00, & 0.0512238276497D+00, & 0.2385615300181D+00, & 0.5228769399639D+00, & 0.2385615300181D+00 /) elseif ( strength == 9 ) then x(1:n) = (/ & 0.0451890097844D+00, & 0.0451890097844D+00, & 0.9096219804312D+00, & 0.7475124727339D+00, & 0.2220631655373D+00, & 0.7475124727339D+00, & 0.2220631655373D+00, & 0.0304243617288D+00, & 0.0304243617288D+00, & 0.1369912012649D+00, & 0.6447187277637D+00, & 0.1369912012649D+00, & 0.2182900709714D+00, & 0.2182900709714D+00, & 0.6447187277637D+00, & 0.0369603304334D+00, & 0.4815198347833D+00, & 0.4815198347833D+00, & 0.4036039798179D+00, & 0.4036039798179D+00, & 0.1927920403641D+00 /) elseif ( strength == 11 ) then x(1:n) = (/ & 0.0000000000000D+00, & 0.9451704450173D+00, & 0.9289002405719D+00, & 0.0685505797224D+00, & 0.0243268355615D+00, & 0.1279662835335D+00, & 0.0277838749488D+00, & 0.0287083428360D+00, & 0.7498347588656D+00, & 0.7228007909707D+00, & 0.2497602062386D+00, & 0.0865562992839D+00, & 0.8325513856998D+00, & 0.3061619157672D+00, & 0.0303526617491D+00, & 0.4868610595047D+00, & 0.6657904293017D+00, & 0.1765456154221D+00, & 0.0293121007360D+00, & 0.5295657488667D+00, & 0.1444673824391D+00, & 0.3299740111411D+00, & 0.5361815729052D+00, & 0.5511507516862D+00, & 0.1437790861923D+00, & 0.3348066587327D+00, & 0.1529619437161D+00, & 0.3430183498147D+00 /) elseif ( strength == 13 ) then x(1:n) = (/ & 0.0242935351590D+00, & 0.0265193427722D+00, & 0.9492126023551D+00, & 0.0033775763749D+00, & 0.4757672298101D+00, & 0.5190783193471D+00, & 0.8616839745321D+00, & 0.1249209759926D+00, & 0.0138565453861D+00, & 0.0211887064222D+00, & 0.8432296787219D+00, & 0.1354231797865D+00, & 0.3088853510679D+00, & 0.6685057595169D+00, & 0.0226545012557D+00, & 0.2808515408772D+00, & 0.6922446749051D+00, & 0.0268617447119D+00, & 0.1141778485470D+00, & 0.7974807922061D+00, & 0.0892807293894D+00, & 0.1052487892455D+00, & 0.6663022280740D+00, & 0.2307803737547D+00, & 0.1705059157540D+00, & 0.5086593973043D+00, & 0.3141823862281D+00, & 0.4617460817864D+00, & 0.0693087496081D+00, & 0.4651955259268D+00, & 0.2578625857893D+00, & 0.6112627766779D+00, & 0.1305182135934D+00, & 0.4281437991828D+00, & 0.3356995783730D+00, & 0.2305424298836D+00 /) elseif ( strength == 14 ) then x(1:n) = (/ & 0.0000000000000D+00, & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0573330873026D+00, & 0.0573330873026D+00, & 0.9275286857160D+00, & 0.0151382269814D+00, & 0.9275286857160D+00, & 0.0151382269814D+00, & 0.8159625040711D+00, & 0.8159625040711D+00, & 0.1659719969565D+00, & 0.0180654989724D+00, & 0.1659719969565D+00, & 0.0180654989724D+00, & 0.3165475556378D+00, & 0.6647637544849D+00, & 0.0186886898773D+00, & 0.0186886898773D+00, & 0.3165475556378D+00, & 0.6647637544849D+00, & 0.0192662192492D+00, & 0.4903668903754D+00, & 0.4903668903754D+00, & 0.0875134669581D+00, & 0.0875134669581D+00, & 0.8249730660837D+00, & 0.0935526036219D+00, & 0.0935526036219D+00, & 0.2079865423167D+00, & 0.6984608540613D+00, & 0.6984608540613D+00, & 0.2079865423167D+00, & 0.0974892983467D+00, & 0.3645018421383D+00, & 0.5380088595149D+00, & 0.5380088595149D+00, & 0.3645018421383D+00, & 0.0974892983467D+00, & 0.2217145894873D+00, & 0.5565708210253D+00, & 0.2217145894873D+00, & 0.3860471669296D+00, & 0.2279056661408D+00, & 0.3860471669296D+00 /) elseif ( strength == 16 ) then x(1:n) = (/ & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0000000000000D+00, & 0.9398863583577D+00, & 0.0543806683058D+00, & 0.0093940049164D+00, & 0.0164345086362D+00, & 0.9469487269862D+00, & 0.0426604005768D+00, & 0.0122269495439D+00, & 0.8673696521047D+00, & 0.8456744021389D+00, & 0.1395759632103D+00, & 0.1317821743231D+00, & 0.0157955126300D+00, & 0.7365462884436D+00, & 0.0139688430330D+00, & 0.2547895186039D+00, & 0.7316386522555D+00, & 0.0157253728951D+00, & 0.2662302843647D+00, & 0.8673504065214D+00, & 0.0741493666957D+00, & 0.0159285948360D+00, & 0.0156061028068D+00, & 0.5910094817484D+00, & 0.4034771496889D+00, & 0.5694745628526D+00, & 0.0678493700650D+00, & 0.4265968590272D+00, & 0.0670982507890D+00, & 0.7528310231480D+00, & 0.7753727783557D+00, & 0.1689073157787D+00, & 0.1687335832919D+00, & 0.0821244708436D+00, & 0.6288705363345D+00, & 0.0811413015266D+00, & 0.2969112065080D+00, & 0.0767542314171D+00, & 0.6223022333845D+00, & 0.3103786288051D+00, & 0.0819218215187D+00, & 0.4717022665013D+00, & 0.4546603415250D+00, & 0.1701091339237D+00, & 0.6406004329487D+00, & 0.1912267583717D+00, & 0.1885315767070D+00, & 0.4772929957691D+00, & 0.3126974621760D+00, & 0.4961225945946D+00, & 0.1928805312867D+00, & 0.3360041453816D+00, & 0.3337280550848D+00 /) elseif ( strength == 18 ) then x(1:n) = (/ & 0.0116731059668D+00, & 0.9810030858388D+00, & 0.0106966317092D+00, & 0.9382476983551D+00, & 0.0126627518417D+00, & 0.0598109409984D+00, & 0.0137363297927D+00, & 0.9229527959405D+00, & 0.0633107354993D+00, & 0.0117265100335D+00, & 0.1554720587323D+00, & 0.8343293888982D+00, & 0.8501638031957D+00, & 0.0128816350522D+00, & 0.1510801608959D+00, & 0.0101917879217D+00, & 0.2813372399303D+00, & 0.7124374628501D+00, & 0.2763025250863D+00, & 0.0109658368561D+00, & 0.4289110517884D+00, & 0.4215420555115D+00, & 0.5711258590444D+00, & 0.5826868270511D+00, & 0.0130567806713D+00, & 0.0130760400964D+00, & 0.7263437062407D+00, & 0.0687230068637D+00, & 0.8652302101529D+00, & 0.0648599071037D+00, & 0.1483494943362D+00, & 0.0624359898396D+00, & 0.7871369011735D+00, & 0.0519104921610D+00, & 0.1543129927444D+00, & 0.2617842745603D+00, & 0.7667257872813D+00, & 0.2582103676627D+00, & 0.0679065925147D+00, & 0.5293578274804D+00, & 0.0666036150484D+00, & 0.0585675461899D+00, & 0.0644535360411D+00, & 0.6748138429151D+00, & 0.3914602310369D+00, & 0.6487701492307D+00, & 0.3946498220408D+00, & 0.5390137151933D+00, & 0.1627895082785D+00, & 0.6812436322641D+00, & 0.1542832878020D+00, & 0.2522727750445D+00, & 0.2547981532407D+00, & 0.1485580549194D+00, & 0.2930239606436D+00, & 0.2808991272310D+00, & 0.4820989592971D+00, & 0.5641878245444D+00, & 0.1307699644344D+00, & 0.1479692221948D+00, & 0.5638684222946D+00, & 0.4361157428790D+00, & 0.3603263935285D+00, & 0.4224188334674D+00, & 0.3719001833052D+00, & 0.2413645006928D+00 /) elseif ( strength == 20 ) then x(1:n) = (/ & 0.0089411337112D+00, & 0.9792622629807D+00, & 0.0105475382112D+00, & 0.0023777061947D+00, & 0.0630425115795D+00, & 0.9308422496730D+00, & 0.0629076555490D+00, & 0.9315962246381D+00, & 0.0061951689415D+00, & 0.0287125819237D+00, & 0.9293844478305D+00, & 0.0375457566621D+00, & 0.0086895739064D+00, & 0.1547597053965D+00, & 0.8331025294185D+00, & 0.8374231073526D+00, & 0.1559362505234D+00, & 0.0098599642095D+00, & 0.4055873733289D+00, & 0.5964727898618D+00, & 0.0080747800416D+00, & 0.0075073977721D+00, & 0.3936764519237D+00, & 0.5846530726212D+00, & 0.4870804112120D+00, & 0.2683512811785D+00, & 0.7223956288748D+00, & 0.2716826742357D+00, & 0.0112580842046D+00, & 0.0115034734370D+00, & 0.7140525900564D+00, & 0.4902871053112D+00, & 0.0201423425209D+00, & 0.0361107464859D+00, & 0.8607998819851D+00, & 0.1005891526001D+00, & 0.0918740717058D+00, & 0.8604888296191D+00, & 0.0439842178673D+00, & 0.2011017606735D+00, & 0.7449993726263D+00, & 0.0532186641310D+00, & 0.7453984647401D+00, & 0.1957289932876D+00, & 0.1092532057988D+00, & 0.0567625702001D+00, & 0.0483837933475D+00, & 0.1080612809760D+00, & 0.6185605900991D+00, & 0.7721296013497D+00, & 0.6115734801133D+00, & 0.3381326103376D+00, & 0.1173084128254D+00, & 0.2674551260596D+00, & 0.6542100160026D+00, & 0.0538297481158D+00, & 0.1848840324117D+00, & 0.3376267104744D+00, & 0.6067102034499D+00, & 0.4612614085496D+00, & 0.1525465365671D+00, & 0.0700582543543D+00, & 0.4704201379032D+00, & 0.1216461693746D+00, & 0.6371404052702D+00, & 0.2379904515119D+00, & 0.1483929857177D+00, & 0.3598069571550D+00, & 0.4941441055095D+00, & 0.1440630687981D+00, & 0.5019764440004D+00, & 0.3555423834298D+00, & 0.2443439540771D+00, & 0.2437064989342D+00, & 0.5122200807321D+00, & 0.2526038315178D+00, & 0.3759895652851D+00, & 0.3729077987144D+00 /) elseif ( strength == 21 ) then x(1:n) = (/ & 0.0035524391922D+00, & 0.0035524391922D+00, & 0.9928951216156D+00, & 0.9553548273730D+00, & 0.0358552797177D+00, & 0.9553548273730D+00, & 0.0087898929093D+00, & 0.0087898929093D+00, & 0.0358552797177D+00, & 0.8865264879047D+00, & 0.8865264879047D+00, & 0.0052405375935D+00, & 0.0052405375935D+00, & 0.1082329745017D+00, & 0.1082329745017D+00, & 0.0466397432150D+00, & 0.0466397432150D+00, & 0.9067205135700D+00, & 0.2075720456946D+00, & 0.2075720456946D+00, & 0.7841520301770D+00, & 0.0082759241284D+00, & 0.0082759241284D+00, & 0.7841520301770D+00, & 0.0858119489725D+00, & 0.8827043562574D+00, & 0.0314836947701D+00, & 0.0858119489725D+00, & 0.8827043562574D+00, & 0.0314836947701D+00, & 0.6688778233826D+00, & 0.0095150760625D+00, & 0.0095150760625D+00, & 0.6688778233826D+00, & 0.3216071005550D+00, & 0.3216071005550D+00, & 0.4379999543113D+00, & 0.0099859785681D+00, & 0.4379999543113D+00, & 0.0099859785681D+00, & 0.5520140671206D+00, & 0.5520140671206D+00, & 0.7974931072148D+00, & 0.0405093994119D+00, & 0.0405093994119D+00, & 0.1619974933734D+00, & 0.7974931072148D+00, & 0.1619974933734D+00, & 0.3864215551955D+00, & 0.3864215551955D+00, & 0.2271568896090D+00, & 0.8090129379329D+00, & 0.0954935310336D+00, & 0.0954935310336D+00, & 0.2745425238718D+00, & 0.0479840480721D+00, & 0.6774734280561D+00, & 0.6774734280561D+00, & 0.2745425238718D+00, & 0.0479840480721D+00, & 0.4053472446667D+00, & 0.0516677930989D+00, & 0.4053472446667D+00, & 0.5429849622344D+00, & 0.0516677930989D+00, & 0.5429849622344D+00, & 0.1877738615539D+00, & 0.7054113116872D+00, & 0.7054113116872D+00, & 0.1068148267588D+00, & 0.1877738615539D+00, & 0.1068148267588D+00, & 0.1195059712009D+00, & 0.1195059712009D+00, & 0.5747817297348D+00, & 0.5747817297348D+00, & 0.3057122990643D+00, & 0.3057122990643D+00, & 0.5981245743363D+00, & 0.2009377128319D+00, & 0.2009377128319D+00, & 0.2160775200005D+00, & 0.3121360256673D+00, & 0.2160775200005D+00, & 0.3121360256673D+00, & 0.4717864543321D+00, & 0.4717864543321D+00, & 0.4376579903849D+00, & 0.4376579903849D+00, & 0.1246840192303D+00, & 0.3333333333333D+00 /) elseif ( strength == 23 ) then x(1:n) = (/ & 0.0087809303836D+00, & 0.9903675314220D+00, & 0.0027029276450D+00, & 0.0335909214524D+00, & 0.0091675068606D+00, & 0.9675568182558D+00, & 0.0084737200688D+00, & 0.0078781948792D+00, & 0.0676785477700D+00, & 0.9470266955047D+00, & 0.0442974755680D+00, & 0.9144243214882D+00, & 0.0081735424459D+00, & 0.2497452292741D+00, & 0.3833232646055D+00, & 0.8876850353557D+00, & 0.1035329228297D+00, & 0.0077255923618D+00, & 0.1403192425107D+00, & 0.8104591009652D+00, & 0.1809643003717D+00, & 0.8330767948684D+00, & 0.0083010907126D+00, & 0.0348407706147D+00, & 0.2740287679608D+00, & 0.7173982224778D+00, & 0.2394976858234D+00, & 0.0081859185845D+00, & 0.0068836152075D+00, & 0.4843741485699D+00, & 0.4960767772741D+00, & 0.6112936776245D+00, & 0.3804323980345D+00, & 0.7303890713524D+00, & 0.0083987168639D+00, & 0.6128525675612D+00, & 0.0075475961037D+00, & 0.0079525316513D+00, & 0.3559774870460D+00, & 0.9110236977966D+00, & 0.0437233605166D+00, & 0.0388480061835D+00, & 0.0967032117936D+00, & 0.0873226911312D+00, & 0.0421445202084D+00, & 0.8485617974961D+00, & 0.8477921333864D+00, & 0.1067435889398D+00, & 0.1833966521991D+00, & 0.0416340541167D+00, & 0.7611632251560D+00, & 0.1941599254144D+00, & 0.7579378747173D+00, & 0.0439826512395D+00, & 0.0369760535918D+00, & 0.5363187134342D+00, & 0.1001256948921D+00, & 0.7912266693524D+00, & 0.0379866714177D+00, & 0.4157414028965D+00, & 0.6507106491463D+00, & 0.0420141133438D+00, & 0.0425548444254D+00, & 0.2920627107240D+00, & 0.5389729538180D+00, & 0.4193031828489D+00, & 0.6549472009700D+00, & 0.3007352790917D+00, & 0.3752400771585D+00, & 0.3453980282786D+00, & 0.0994532168761D+00, & 0.1598309359585D+00, & 0.1797326661667D+00, & 0.7124584461943D+00, & 0.1066065678636D+00, & 0.7001701904096D+00, & 0.0993303629801D+00, & 0.6065648052521D+00, & 0.1023223542704D+00, & 0.2533382324938D+00, & 0.6166226715217D+00, & 0.2769500693109D+00, & 0.0904184571873D+00, & 0.4981522767248D+00, & 0.0928231860168D+00, & 0.3738418699229D+00, & 0.2521678840407D+00, & 0.5087500218708D+00, & 0.3905579116731D+00, & 0.1706141469096D+00, & 0.5266737761312D+00, & 0.3487581527629D+00, & 0.2588053596017D+00, & 0.1696614558053D+00, & 0.3013521806875D+00, & 0.2580202409759D+00, & 0.4584740860198D+00, & 0.1848898683498D+00, & 0.6130740338465D+00, & 0.1921611750994D+00, & 0.4180541160599D+00, & 0.1650612642036D+00, & 0.5159205739625D+00, & 0.2982718935750D+00, & 0.4098894602340D+00 /) elseif ( strength == 25 ) then x(1:n) = (/ & 0.0082881595033D+00, & 0.4618422030241D+00, & 0.0071066441239D+00, & 0.9847613141699D+00, & 0.5374447869049D+00, & 0.0000000000000D+00, & 0.4914131929361D+00, & 0.0070345937020D+00, & 0.9564734714228D+00, & 0.0370198792045D+00, & 0.1024124542747D+00, & 0.5928065811509D+00, & 0.0050948422371D+00, & 0.0081562023689D+00, & 0.0424936107568D+00, & 0.9495543500844D+00, & 0.8932787471239D+00, & 0.0069317612927D+00, & 0.9035839030665D+00, & 0.0905665738209D+00, & 0.0083929332787D+00, & 0.6261245686071D+00, & 0.0062801592979D+00, & 0.8272539257367D+00, & 0.0062005875353D+00, & 0.1676900311185D+00, & 0.7199353069567D+00, & 0.2749740090237D+00, & 0.0079257582005D+00, & 0.0069981220752D+00, & 0.8125248773263D+00, & 0.0073536969970D+00, & 0.7283665935411D+00, & 0.1800642304565D+00, & 0.2658102467762D+00, & 0.0070892364520D+00, & 0.3774054302043D+00, & 0.0369649608668D+00, & 0.9203194109805D+00, & 0.0425477806431D+00, & 0.6191278394983D+00, & 0.3762697209178D+00, & 0.0956111149690D+00, & 0.0302473410377D+00, & 0.8739905691754D+00, & 0.8604133734958D+00, & 0.0347307852352D+00, & 0.1043606608343D+00, & 0.7797622824754D+00, & 0.0185865164256D+00, & 0.0324585286618D+00, & 0.8371293901157D+00, & 0.0836602075315D+00, & 0.0784070242501D+00, & 0.4929238648458D+00, & 0.1870637584073D+00, & 0.4892636967025D+00, & 0.0401982618372D+00, & 0.7894259278865D+00, & 0.1686260456429D+00, & 0.3750901913174D+00, & 0.0356362876880D+00, & 0.5887548164804D+00, & 0.0373308082182D+00, & 0.2820769993374D+00, & 0.6819277603320D+00, & 0.0374938324382D+00, & 0.6984079204127D+00, & 0.2654390894079D+00, & 0.1429848440800D+00, & 0.7623554007647D+00, & 0.0934222022749D+00, & 0.5759004479923D+00, & 0.3822427332525D+00, & 0.0411414081675D+00, & 0.0802462538379D+00, & 0.7625229819410D+00, & 0.1524941445131D+00, & 0.0622159195833D+00, & 0.1109539036076D+00, & 0.4575627212057D+00, & 0.4322865136374D+00, & 0.5865002850241D+00, & 0.0869359250818D+00, & 0.0929594906936D+00, & 0.6661932141454D+00, & 0.4780306362227D+00, & 0.4372215294577D+00, & 0.6779224504669D+00, & 0.2423431255660D+00, & 0.2288925420305D+00, & 0.3315065049959D+00, & 0.3424200526607D+00, & 0.0862630046475D+00, & 0.5113188946635D+00, & 0.1538977841001D+00, & 0.6779951348472D+00, & 0.1664600469411D+00, & 0.0950910318888D+00, & 0.3436048136712D+00, & 0.5560417025366D+00, & 0.1452404029513D+00, & 0.1619685156238D+00, & 0.5800164844262D+00, & 0.2450201223288D+00, & 0.2557621891794D+00, & 0.2205239985511D+00, & 0.4940183111285D+00, & 0.2531570689798D+00, & 0.5846891116357D+00, & 0.1660333602278D+00, & 0.2505426292461D+00, & 0.3519336802182D+00, & 0.3502668835419D+00, & 0.4400892485512D+00, & 0.4680855471546D+00, & 0.1770237763947D+00, & 0.3900920779501D+00, & 0.2805847774120D+00, & 0.3361523347440D+00 /) end if return end subroutine twb_rule_y ( strength, y ) !*****************************************************************************80 ! !! twb_rule_y returns the y abscissas of a TWB rule of given strength. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 April 2019 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Mark Taylor, Beth Wingate, Len Bos, ! Several new quadrature formulas for polynomial integration in the triangle, ! http://arxiv.org/abs/math/0501496v2, ! 08 February 2007. ! ! Parameters: ! ! Input, integer STRENGTH, the desired strength. ! 1 <= STRENGTH. ! ! Output, real ( kind = rk ) Y(N), the y abscissas of the rule, if it exists. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer strength real ( kind = rk ) y(*) call twb_rule_n ( strength, n ) if ( n < 1 ) then return end if if ( strength == 1 ) then y(1:n) = (/ & 0.3333333333333D+00 /) else if ( strength == 2 ) then y(1:n) = (/ & 0.6666666666667D+00, & 0.1666666666667D+00, & 0.1666666666667D+00 /) else if ( strength == 4 ) then y(1:n) = (/ & 0.0915762135098D+00, & 0.0915762135098D+00, & 0.8168475729805D+00, & 0.4459484909160D+00, & 0.1081030181681D+00, & 0.4459484909160D+00 /) else if ( strength == 5 ) then y(1:n) = (/ & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0000000000000D+00, & 0.6728199218710D+00, & 0.2673288599482D+00, & 0.6716530111494D+00, & 0.0649251690029D+00, & 0.2693789366453D+00, & 0.0654054874919D+00, & 0.3386799893027D+00 /) else if ( strength == 7 ) then y(1:n) = (/ & 0.0000000000000D+00, & 0.0000000000000D+00, & 1.0000000000000D+00, & 0.0421382841642D+00, & 0.7839656651012D+00, & 0.0421382841642D+00, & 0.1738960507345D+00, & 0.1738960507345D+00, & 0.7839656651012D+00, & 0.4743880861752D+00, & 0.0512238276497D+00, & 0.4743880861752D+00, & 0.5228769399639D+00, & 0.2385615300181D+00, & 0.2385615300181D+00 /) else if ( strength == 9 ) then y(1:n) = (/ & 0.0451890097844D+00, & 0.9096219804312D+00, & 0.0451890097844D+00, & 0.0304243617288D+00, & 0.0304243617288D+00, & 0.2220631655373D+00, & 0.7475124727339D+00, & 0.7475124727339D+00, & 0.2220631655373D+00, & 0.2182900709714D+00, & 0.2182900709714D+00, & 0.6447187277637D+00, & 0.6447187277637D+00, & 0.1369912012649D+00, & 0.1369912012649D+00, & 0.4815198347833D+00, & 0.0369603304334D+00, & 0.4815198347833D+00, & 0.1927920403641D+00, & 0.4036039798179D+00, & 0.4036039798179D+00 /) else if ( strength == 11 ) then y(1:n) = (/ & 0.9451704450174D+00, & 0.0000000000000D+00, & 0.0685505797224D+00, & 0.9289002405717D+00, & 0.0243268355616D+00, & 0.0277838749488D+00, & 0.1279662835337D+00, & 0.7498347588657D+00, & 0.0287083428360D+00, & 0.2497602062385D+00, & 0.7228007909707D+00, & 0.8325513856997D+00, & 0.0865562992839D+00, & 0.0303526617491D+00, & 0.3061619157675D+00, & 0.4868610595047D+00, & 0.1765456154219D+00, & 0.6657904293016D+00, & 0.5295657488669D+00, & 0.0293121007360D+00, & 0.1444673824391D+00, & 0.5361815729050D+00, & 0.3299740111409D+00, & 0.1437790861923D+00, & 0.5511507516862D+00, & 0.1529619437161D+00, & 0.3348066587327D+00, & 0.3430183498147D+00 /) else if ( strength == 13 ) then y(1:n) = (/ & 0.9493059293846D+00, & 0.0242695130640D+00, & 0.0265067966437D+00, & 0.4767316412363D+00, & 0.5198921829102D+00, & 0.0055912706202D+00, & 0.0133996048618D+00, & 0.8613054321334D+00, & 0.1247733717358D+00, & 0.8438438351223D+00, & 0.1354563645830D+00, & 0.0213482820656D+00, & 0.0221919663014D+00, & 0.3089012879389D+00, & 0.6691709943321D+00, & 0.6924718155106D+00, & 0.0268723345026D+00, & 0.2810093973222D+00, & 0.7973581413586D+00, & 0.0879806508791D+00, & 0.1145020561128D+00, & 0.6686904119922D+00, & 0.2275051631832D+00, & 0.1054572561221D+00, & 0.5174064398658D+00, & 0.3170523855209D+00, & 0.1810706361659D+00, & 0.4678594539804D+00, & 0.4622856042085D+00, & 0.0724357805669D+00, & 0.6131395039177D+00, & 0.1300360834609D+00, & 0.2581713828884D+00, & 0.2362005969817D+00, & 0.4311026308588D+00, & 0.3456013949376D+00 /) else if ( strength == 14 ) then y(1:n) = (/ & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0000000000000D+00, & 0.0151382269814D+00, & 0.9275286857160D+00, & 0.0573330873026D+00, & 0.0573330873026D+00, & 0.0151382269814D+00, & 0.9275286857160D+00, & 0.1659719969565D+00, & 0.0180654989724D+00, & 0.8159625040711D+00, & 0.8159625040711D+00, & 0.0180654989724D+00, & 0.1659719969565D+00, & 0.0186886898773D+00, & 0.0186886898773D+00, & 0.6647637544849D+00, & 0.3165475556378D+00, & 0.6647637544849D+00, & 0.3165475556378D+00, & 0.4903668903754D+00, & 0.0192662192492D+00, & 0.4903668903754D+00, & 0.8249730660837D+00, & 0.0875134669581D+00, & 0.0875134669581D+00, & 0.2079865423167D+00, & 0.6984608540613D+00, & 0.0935526036219D+00, & 0.0935526036219D+00, & 0.2079865423167D+00, & 0.6984608540613D+00, & 0.5380088595149D+00, & 0.0974892983467D+00, & 0.0974892983467D+00, & 0.3645018421383D+00, & 0.5380088595149D+00, & 0.3645018421383D+00, & 0.5565708210253D+00, & 0.2217145894873D+00, & 0.2217145894873D+00, & 0.2279056661408D+00, & 0.3860471669296D+00, & 0.3860471669296D+00 /) else if ( strength == 16 ) then y(1:n) = (/ & 0.0000000000000D+00, & 1.0000000000000D+00, & 0.0000000000000D+00, & 0.0049848744634D+00, & 0.9386405618617D+00, & 0.0526424462697D+00, & 0.9469035517351D+00, & 0.0363373677167D+00, & 0.0151224541799D+00, & 0.8693773510664D+00, & 0.1204917285774D+00, & 0.0157763967870D+00, & 0.8448120870375D+00, & 0.0135009605584D+00, & 0.1455274938536D+00, & 0.0155697540908D+00, & 0.7379836894450D+00, & 0.7297615689771D+00, & 0.2543076683315D+00, & 0.2696239795791D+00, & 0.0144783956308D+00, & 0.0591679410400D+00, & 0.8634782575061D+00, & 0.4191238955238D+00, & 0.5809222921146D+00, & 0.0159251452651D+00, & 0.5806700368104D+00, & 0.4149495146302D+00, & 0.0761218678591D+00, & 0.0157509692312D+00, & 0.7741898312421D+00, & 0.0819119495639D+00, & 0.1577128457292D+00, & 0.7503943099742D+00, & 0.0708311507268D+00, & 0.1762996626771D+00, & 0.0807744953317D+00, & 0.3054373589776D+00, & 0.6227485988871D+00, & 0.6247247149546D+00, & 0.3011485821166D+00, & 0.0779098365079D+00, & 0.4603633038351D+00, & 0.0821554006797D+00, & 0.4637565033890D+00, & 0.6422277808188D+00, & 0.1898293537256D+00, & 0.1739955685343D+00, & 0.4798914070406D+00, & 0.3348356598119D+00, & 0.4957972197259D+00, & 0.1927553668904D+00, & 0.3161015807261D+00, & 0.1894892801290D+00, & 0.3343571021811D+00 /) else if ( strength == 18 ) then y(1:n) = (/ & 0.9812565951289D+00, & 0.0071462504863D+00, & 0.0115153933376D+00, & 0.0495570591341D+00, & 0.9370123620615D+00, & 0.0121364578922D+00, & 0.0612783625597D+00, & 0.0141128270602D+00, & 0.9220197291727D+00, & 0.1500520475229D+00, & 0.8325147121589D+00, & 0.0125228158759D+00, & 0.1371997508736D+00, & 0.8477627063479D+00, & 0.0136526924039D+00, & 0.5770438618345D+00, & 0.7066853759623D+00, & 0.0124569780990D+00, & 0.0121741311386D+00, & 0.4194306712466D+00, & 0.5599616067469D+00, & 0.0116475994785D+00, & 0.0118218313989D+00, & 0.4057889581177D+00, & 0.2725023750868D+00, & 0.7224712523233D+00, & 0.2602984019251D+00, & 0.0631417277210D+00, & 0.0720611837338D+00, & 0.8590433543910D+00, & 0.7888788352240D+00, & 0.1493935499354D+00, & 0.0656382042757D+00, & 0.5255635695605D+00, & 0.0716383926917D+00, & 0.0621479485288D+00, & 0.1658211554831D+00, & 0.6800119766139D+00, & 0.7571515437782D+00, & 0.4121503841107D+00, & 0.2612513087886D+00, & 0.3902236114535D+00, & 0.6373626559761D+00, & 0.0637583342061D+00, & 0.5503238090563D+00, & 0.2836728360263D+00, & 0.0605175522554D+00, & 0.0611990176936D+00, & 0.6861322141035D+00, & 0.1567968345899D+00, & 0.1667512624020D+00, & 0.2504803933395D+00, & 0.4994090649043D+00, & 0.5756023096087D+00, & 0.5656897354162D+00, & 0.1437921574248D+00, & 0.2518557535865D+00, & 0.1462966743153D+00, & 0.4489577586117D+00, & 0.3001174386829D+00, & 0.2813772089298D+00, & 0.4252053446420D+00, & 0.2599190004889D+00, & 0.1453238443303D+00, & 0.3780122703567D+00, & 0.3847563284940D+00 /) else if ( strength == 20 ) then y(1:n) = (/ & 0.0086983293702D+00, & 0.0102644133744D+00, & 0.9785514202515D+00, & 0.0636551098604D+00, & 0.0041506347509D+00, & 0.0048053482263D+00, & 0.9316790069481D+00, & 0.0626264881801D+00, & 0.9293587058564D+00, & 0.0310202122997D+00, & 0.0342152968219D+00, & 0.9257868884669D+00, & 0.1584971251510D+00, & 0.8363606657688D+00, & 0.0089257244824D+00, & 0.1529167304078D+00, & 0.0094966240058D+00, & 0.8342211493596D+00, & 0.0074389302008D+00, & 0.3956330809311D+00, & 0.4031319425903D+00, & 0.5851609594681D+00, & 0.5974896592899D+00, & 0.0087250464968D+00, & 0.0202129229912D+00, & 0.7202340088668D+00, & 0.2662399366456D+00, & 0.0112882698808D+00, & 0.7169695963325D+00, & 0.2740067110166D+00, & 0.0113511560497D+00, & 0.4936491841468D+00, & 0.4832573459601D+00, & 0.0935679501582D+00, & 0.0397379067075D+00, & 0.8586343419352D+00, & 0.0395513001973D+00, & 0.0966224057079D+00, & 0.8561886349107D+00, & 0.7449115835626D+00, & 0.0536865638166D+00, & 0.1963754275935D+00, & 0.1982065805550D+00, & 0.0555713833156D+00, & 0.6100036182413D+00, & 0.7409121894959D+00, & 0.6075135660978D+00, & 0.1122081510437D+00, & 0.2698753703035D+00, & 0.1114117395333D+00, & 0.3389367677931D+00, & 0.0494693938787D+00, & 0.7696451309795D+00, & 0.1115718808154D+00, & 0.1906548314700D+00, & 0.3358616826849D+00, & 0.1551831523851D+00, & 0.6081402596294D+00, & 0.0542632795598D+00, & 0.0688176670722D+00, & 0.6510240845749D+00, & 0.4661904392742D+00, & 0.4634826455353D+00, & 0.2381494875516D+00, & 0.1238399384513D+00, & 0.6370216452326D+00, & 0.4894188577780D+00, & 0.1452880866253D+00, & 0.3610216383818D+00, & 0.3513508341887D+00, & 0.1435491663293D+00, & 0.5016491599502D+00, & 0.2406052129104D+00, & 0.5109017277055D+00, & 0.2452737973543D+00, & 0.3700319555094D+00, & 0.2505406611631D+00, & 0.3753750277549D+00 /) else if ( strength == 21 ) then y(1:n) = (/ & 0.0035524391922D+00, & 0.9928951216156D+00, & 0.0035524391922D+00, & 0.0087898929093D+00, & 0.0087898929093D+00, & 0.0358552797177D+00, & 0.0358552797177D+00, & 0.9553548273730D+00, & 0.9553548273730D+00, & 0.1082329745017D+00, & 0.0052405375935D+00, & 0.1082329745017D+00, & 0.8865264879047D+00, & 0.8865264879047D+00, & 0.0052405375935D+00, & 0.9067205135700D+00, & 0.0466397432150D+00, & 0.0466397432150D+00, & 0.0082759241284D+00, & 0.7841520301770D+00, & 0.2075720456946D+00, & 0.7841520301770D+00, & 0.2075720456946D+00, & 0.0082759241284D+00, & 0.0314836947701D+00, & 0.0314836947701D+00, & 0.0858119489725D+00, & 0.8827043562574D+00, & 0.0858119489725D+00, & 0.8827043562574D+00, & 0.0095150760625D+00, & 0.3216071005550D+00, & 0.6688778233826D+00, & 0.3216071005550D+00, & 0.6688778233826D+00, & 0.0095150760625D+00, & 0.0099859785681D+00, & 0.5520140671206D+00, & 0.5520140671206D+00, & 0.4379999543113D+00, & 0.4379999543113D+00, & 0.0099859785681D+00, & 0.0405093994119D+00, & 0.1619974933734D+00, & 0.7974931072148D+00, & 0.7974931072148D+00, & 0.1619974933734D+00, & 0.0405093994119D+00, & 0.3864215551955D+00, & 0.2271568896090D+00, & 0.3864215551955D+00, & 0.0954935310336D+00, & 0.8090129379329D+00, & 0.0954935310336D+00, & 0.0479840480721D+00, & 0.6774734280561D+00, & 0.0479840480721D+00, & 0.2745425238718D+00, & 0.6774734280561D+00, & 0.2745425238718D+00, & 0.5429849622344D+00, & 0.4053472446667D+00, & 0.0516677930989D+00, & 0.0516677930989D+00, & 0.5429849622344D+00, & 0.4053472446667D+00, & 0.1068148267588D+00, & 0.1877738615539D+00, & 0.1068148267588D+00, & 0.7054113116872D+00, & 0.7054113116872D+00, & 0.1877738615539D+00, & 0.3057122990643D+00, & 0.5747817297348D+00, & 0.1195059712009D+00, & 0.3057122990643D+00, & 0.5747817297348D+00, & 0.1195059712009D+00, & 0.2009377128319D+00, & 0.5981245743363D+00, & 0.2009377128319D+00, & 0.3121360256673D+00, & 0.2160775200005D+00, & 0.4717864543321D+00, & 0.4717864543321D+00, & 0.3121360256673D+00, & 0.2160775200005D+00, & 0.4376579903849D+00, & 0.1246840192303D+00, & 0.4376579903849D+00, & 0.3333333333333D+00 /) else if ( strength == 23 ) then y(1:n) = (/ & 0.9903676436772D+00, & 0.0087809216232D+00, & 0.0335914404439D+00, & 0.0027028946710D+00, & 0.0091676353051D+00, & 0.0084737176656D+00, & 0.9675569435345D+00, & 0.0676784943862D+00, & 0.0078781659291D+00, & 0.0442974541187D+00, & 0.9470266676487D+00, & 0.0081735455132D+00, & 0.9144244234031D+00, & 0.3833232434720D+00, & 0.2497451268005D+00, & 0.1035328809446D+00, & 0.8876849931840D+00, & 0.1403190991974D+00, & 0.0077255934624D+00, & 0.1809642523926D+00, & 0.8104590515334D+00, & 0.0083010939677D+00, & 0.8330768545392D+00, & 0.0348406969482D+00, & 0.7173981847948D+00, & 0.2740287304386D+00, & 0.0081859182262D+00, & 0.2394975566677D+00, & 0.4843740892687D+00, & 0.0068836232949D+00, & 0.4960767529507D+00, & 0.3804323691239D+00, & 0.6112936466533D+00, & 0.0083987179701D+00, & 0.7303890895407D+00, & 0.0075475979695D+00, & 0.6128525484582D+00, & 0.3559773826721D+00, & 0.0079525358502D+00, & 0.0437233665345D+00, & 0.9110236807446D+00, & 0.0967030908282D+00, & 0.0388479942386D+00, & 0.0873226620391D+00, & 0.8485617789108D+00, & 0.0421445420915D+00, & 0.1067435942472D+00, & 0.8477921328146D+00, & 0.0416340521608D+00, & 0.1833965196930D+00, & 0.1941599202852D+00, & 0.7611632153938D+00, & 0.0439826608586D+00, & 0.7579378242308D+00, & 0.5363186076436D+00, & 0.0369760780935D+00, & 0.7912267093545D+00, & 0.1001257554673D+00, & 0.4157413128558D+00, & 0.0379867061535D+00, & 0.0420141226713D+00, & 0.6507105645084D+00, & 0.2920626023484D+00, & 0.0425548546753D+00, & 0.4193031469005D+00, & 0.5389729093610D+00, & 0.3007352636162D+00, & 0.6549471812731D+00, & 0.3453980130752D+00, & 0.3752400695673D+00, & 0.1598308695187D+00, & 0.0994531960132D+00, & 0.7124585430924D+00, & 0.1797327722240D+00, & 0.7001701784175D+00, & 0.1066065855677D+00, & 0.6065647984796D+00, & 0.0993303896769D+00, & 0.2533381579528D+00, & 0.1023223826189D+00, & 0.2769502060575D+00, & 0.6166227900624D+00, & 0.4981522637001D+00, & 0.0904185045149D+00, & 0.3738418516908D+00, & 0.0928232584790D+00, & 0.2521680925697D+00, & 0.3905580544330D+00, & 0.5087501437661D+00, & 0.5266738039554D+00, & 0.1706142257537D+00, & 0.2588055084886D+00, & 0.3487583491703D+00, & 0.3013522183964D+00, & 0.1696615963219D+00, & 0.4584741774478D+00, & 0.2580203819011D+00, & 0.1848898704551D+00, & 0.1921611994069D+00, & 0.6130740398389D+00, & 0.1650613336416D+00, & 0.4180541199244D+00, & 0.2982719005229D+00, & 0.5159205534362D+00, & 0.4098894317792D+00 /) else if ( strength == 25 ) then y(1:n) = (/ & 0.9848202768869D+00, & 0.5381577969759D+00, & 0.0080842361390D+00, & 0.0070015755134D+00, & 0.4625552130951D+00, & 0.4887676880140D+00, & 0.0000000000000D+00, & 0.9574158053697D+00, & 0.0364655449485D+00, & 0.0070908577166D+00, & 0.8936125594937D+00, & 0.0049451705600D+00, & 0.0996676659189D+00, & 0.0415561148784D+00, & 0.9494865260352D+00, & 0.0081794507292D+00, & 0.0053224326262D+00, & 0.9065401020433D+00, & 0.0894771171077D+00, & 0.0070525342005D+00, & 0.6663179931111D+00, & 0.0092197583153D+00, & 0.8335207460527D+00, & 0.1665134939330D+00, & 0.7424693255229D+00, & 0.0065717743528D+00, & 0.0064354534962D+00, & 0.7185296120719D+00, & 0.1766411374714D+00, & 0.2704767254004D+00, & 0.0082299533210D+00, & 0.5934167875453D+00, & 0.2648817553752D+00, & 0.8115848976682D+00, & 0.0068553525429D+00, & 0.3757632659744D+00, & 0.6148573533757D+00, & 0.9210792302893D+00, & 0.0426025082114D+00, & 0.0372689941794D+00, & 0.3724055713809D+00, & 0.0081436422011D+00, & 0.8771098372601D+00, & 0.0943858903393D+00, & 0.0313198990883D+00, & 0.1049019782046D+00, & 0.8609856462886D+00, & 0.0357152881004D+00, & 0.1872318199265D+00, & 0.4834397678794D+00, & 0.7783474916042D+00, & 0.0804060570156D+00, & 0.8421414817051D+00, & 0.0849927089145D+00, & 0.4892855914710D+00, & 0.0345210858281D+00, & 0.0190774755077D+00, & 0.1691143187109D+00, & 0.0412206731484D+00, & 0.7894860640585D+00, & 0.5895318272013D+00, & 0.3681256217699D+00, & 0.0359968962541D+00, & 0.6790704673533D+00, & 0.0373639992361D+00, & 0.2803330345725D+00, & 0.2634016180014D+00, & 0.0364154673322D+00, & 0.6980717436193D+00, & 0.7612254618453D+00, & 0.0943741220275D+00, & 0.1479799836832D+00, & 0.3821329641698D+00, & 0.0426716362301D+00, & 0.5718082874432D+00, & 0.7702204382042D+00, & 0.1559420577362D+00, & 0.0842965421322D+00, & 0.4538181318873D+00, & 0.4586014071171D+00, & 0.4795313560210D+00, & 0.1230591237472D+00, & 0.0834119779793D+00, & 0.6755677013351D+00, & 0.2326500892727D+00, & 0.2448294007406D+00, & 0.0661749044835D+00, & 0.4442145585244D+00, & 0.0929096534577D+00, & 0.0889793655129D+00, & 0.6780053081672D+00, & 0.5847381559741D+00, & 0.5139245722736D+00, & 0.3340976249234D+00, & 0.1380154720554D+00, & 0.6788062619562D+00, & 0.1663358925269D+00, & 0.1582214504849D+00, & 0.5666590332543D+00, & 0.0978960873457D+00, & 0.3468917820947D+00, & 0.3599534491052D+00, & 0.5810131373330D+00, & 0.2560674640672D+00, & 0.5881469552102D+00, & 0.1652244065047D+00, & 0.3496507466106D+00, & 0.2549448448453D+00, & 0.2543369115017D+00, & 0.1666603916479D+00, & 0.2523240191705D+00, & 0.4959007627528D+00, & 0.1805380367800D+00, & 0.4358582329881D+00, & 0.2120576104941D+00, & 0.3552681570774D+00, & 0.4670352922266D+00, & 0.3323152819300D+00, & 0.3898041176680D+00, & 0.2778500044356D+00 /) end if return end