1 // This code has been mechanically translated from the original FORTRAN 2 // code at http://netlib.org/quadpack. 3 4 /** Authors: Lars Tandle Kyllingstad 5 Copyright: Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved. 6 License: Boost License 1.0 7 */ 8 module scid.ports.quadpack.qmomo; 9 10 11 import std.math; 12 import scid.core.fortran; 13 14 15 16 /// 17 void qmomo(Real)(Real alfa, Real beta, Real* ri_, Real* rj_, Real* rg_, 18 Real* rh_, int integr) 19 { 20 //***begin prologue dqmomo 21 //***date written 820101 (yymmdd) 22 //***revision date 830518 (yymmdd) 23 //***category no. h2a2a1,c3a2 24 //***keywords modified chebyshev moments 25 //***author piessens,robert,appl. math. & progr. div. - k.u.leuven 26 // de doncker,elise,appl. math. & progr. div. - k.u.leuven 27 //***purpose this routine computes modified chebsyshev moments. the k-th 28 // modified chebyshev moment is defined as the integral over 29 // (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev 30 // polynomial of degree k. 31 //***description 32 // 33 // modified chebyshev moments 34 // standard fortran subroutine 35 // double precision version 36 // 37 // parameters 38 // alfa - double precision 39 // parameter in the weight function w(x), alfa.gt.(-1) 40 // 41 // beta - double precision 42 // parameter in the weight function w(x), beta.gt.(-1) 43 // 44 // ri - double precision 45 // vector of dimension 25 46 // ri(k) is the integral over (-1,1) of 47 // (1+x)**alfa*t(k-1,x), k = 1, ..., 25. 48 // 49 // rj - double precision 50 // vector of dimension 25 51 // rj(k) is the integral over (-1,1) of 52 // (1-x)**beta*t(k-1,x), k = 1, ..., 25. 53 // 54 // rg - double precision 55 // vector of dimension 25 56 // rg(k) is the integral over (-1,1) of 57 // (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25. 58 // 59 // rh - double precision 60 // vector of dimension 25 61 // rh(k) is the integral over (-1,1) of 62 // (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25. 63 // 64 // integr - integer 65 // input parameter indicating the modified 66 // moments to be computed 67 // integr = 1 compute ri, rj 68 // = 2 compute ri, rj, rg 69 // = 3 compute ri, rj, rh 70 // = 4 compute ri, rj, rg, rh 71 // 72 //***references (none) 73 //***routines called (none) 74 //***end prologue dqmomo 75 // 76 Real alfp1,alfp2,an,anm1,betp1,betp2,ralf,rbet; 77 int i,im1; 78 // 79 auto rg = dimension(rg_, 25); 80 auto rh = dimension(rh_, 25); 81 auto ri = dimension(ri_, 25); 82 auto rj = dimension(rj_, 25); 83 // 84 // 85 //***first executable statement dqmomo 86 alfp1 = alfa+0.1e1; 87 betp1 = beta+0.1e1; 88 alfp2 = alfa+0.2e1; 89 betp2 = beta+0.2e1; 90 ralf = (cast(Real)0.2e1)^^alfp1; 91 rbet = (cast(Real)0.2e1)^^betp1; 92 // 93 // compute ri, rj using a forward recurrence relation. 94 // 95 ri[1] = ralf/alfp1; 96 rj[1] = rbet/betp1; 97 ri[2] = ri[1]*alfa/alfp2; 98 rj[2] = rj[1]*beta/betp2; 99 an = 0.2e1; 100 anm1 = 0.1e1; 101 for (i=3; i<=25; i++) { //do 20 i=3,25 102 ri[i] = -(ralf+an*(an-alfp2)*ri[i-1])/(anm1*(an+alfp1)); 103 rj[i] = -(rbet+an*(an-betp2)*rj[i-1])/(anm1*(an+betp1)); 104 anm1 = an; 105 an = an+0.1e1; 106 l20: ;} 107 if(integr == 1) goto l70; 108 if(integr == 3) goto l40; 109 // 110 // compute rg using a forward recurrence relation. 111 // 112 rg[1] = -ri[1]/alfp1; 113 rg[2] = -(ralf+ralf)/(alfp2*alfp2)-rg[1]; 114 an = 0.2e1; 115 anm1 = 0.1e1; 116 im1 = 2; 117 for (i=3; i<=25; i++) { //do 30 i=3,25 118 rg[i] = -(an*(an-alfp2)*rg[im1]-an*ri[im1]+anm1*ri[i])/ 119 (anm1*(an+alfp1)); 120 anm1 = an; 121 an = an+0.1e1; 122 im1 = i; 123 l30: ;} 124 if(integr == 2) goto l70; 125 // 126 // compute rh using a forward recurrence relation. 127 // 128 l40: rh[1] = -rj[1]/betp1; 129 rh[2] = -(rbet+rbet)/(betp2*betp2)-rh[1]; 130 an = 0.2e1; 131 anm1 = 0.1e1; 132 im1 = 2; 133 for (i=3; i<=25; i++) { //do 50 i=3,25 134 rh[i] = -(an*(an-betp2)*rh[im1]-an*rj[im1]+ 135 anm1*rj[i])/(anm1*(an+betp1)); 136 anm1 = an; 137 an = an+0.1e1; 138 im1 = i; 139 l50: ;} 140 for (i=2; i<=25; i+=2) { //do 60 i=2,25,2 141 rh[i] = -rh[i]; 142 l60: ;} 143 l70: for (i=2; i<=25; i+=2) { //do 80 i=2,25,2 144 rj[i] = -rj[i]; 145 l80: ;} 146 l90: return; 147 } 148 149 unittest 150 { 151 alias qmomo!float fqmomo; 152 alias qmomo!double dqmomo; 153 alias qmomo!real rqmomo; 154 }