summaryrefslogtreecommitdiff
path: root/Utilib/src/ALGPT.f
diff options
context:
space:
mode:
Diffstat (limited to 'Utilib/src/ALGPT.f')
-rw-r--r--Utilib/src/ALGPT.f408
1 files changed, 408 insertions, 0 deletions
diff --git a/Utilib/src/ALGPT.f b/Utilib/src/ALGPT.f
new file mode 100644
index 0000000..ae46666
--- /dev/null
+++ b/Utilib/src/ALGPT.f
@@ -0,0 +1,408 @@
+*DECK ALGPT
+ SUBROUTINE ALGPT(NGPT,XINF,XSUP,ZGKSI,WGKSI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* returns gauss integration points and weights for integration limits
+* specified and to the order specified
+*
+*Copyright:
+* Copyright (C) 1991 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): R. Roy and O. Vagner
+*
+*Reference:
+* integral of F(X)DX , from XINF to XSUP is given by summ from i=1 to
+* NGPT of F(ZGKSI(I))*WGKSI(I)
+* the points ZGKSI and weights WGKSI are generated according
+* to the technique described in Handbook of mathematical functions
+* M. Abramowitz and I. Stegun, Dover Publication Inc. (1972);
+* all points and weights are taken from A.H. Stroud and D. Secrest,
+* gaussian quadrature formulas, Prentice-Hall (1966).
+*
+*Parameters: input
+* NGPT number of gauss points
+* XINF lower integration limit
+* XSUP upper integration limit
+*
+*Parameters: output
+* ZGKSI integration points
+* WGKSI integration weights
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NGPT
+ REAL XINF,XSUP,ZGKSI(NGPT),WGKSI(NGPT)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER ( I2= 1, I3= 2, I4= 3, I5= 5, I6= 7, I7=10, I8=13,
+ > I9=17,I10=21,I11=26,I12=31,I13=37,I14=43,I15=50,
+ > I16=57,I17=65,I18=73,I19=82,I20=91,I24=101,I28=113,
+ > I32=127,I64=143,IFN=174)
+ PARAMETER ( HALF=0.5E0, TWO=2.0E0 )
+ REAL RZKSI(IFN),RWKSI(IFN),FW,FC
+ CHARACTER CERROR*4
+ SAVE RZKSI,RWKSI
+C N=2
+ DATA (RZKSI(I),I=I2,I3-1) / .577350269189626E0/
+ DATA (RWKSI(I),I=I2,I3-1) / 1.000000000000000E0/
+C N=3
+ DATA (RZKSI(I),I=I3,I4-1) / .774596669241483E0/
+ DATA (RWKSI(I),I=I3,I4-1) / .555555555555556E0/
+C N=4
+ DATA (RZKSI(I),I=I4,I5-1) / .339981043584856E0,
+ > .861136311594053E0/
+ DATA (RWKSI(I),I=I4,I5-1) / .652145154862546E0,
+ > .347854845137454E0/
+C N=5
+ DATA (RZKSI(I),I=I5,I6-1) / .538469310105683E0,
+ > .906179845938664E0/
+ DATA (RWKSI(I),I=I5,I6-1) / .478628670499366E0,
+ > .236926885056189E0/
+C N=6
+ DATA (RZKSI(I),I=I6,I7-1) / .238619186083197E0,
+ > .661209386466265E0,
+ > .932469514203152E0/
+ DATA (RWKSI(I),I=I6,I7-1) / .467913934572691E0,
+ > .360761573048139E0,
+ > .171324492379170E0/
+C N=7
+ DATA (RZKSI(I),I=I7,I8-1) / .405845151377397E0,
+ > .741531185599394E0,
+ > .949107912342759E0/
+ DATA (RWKSI(I),I=I7,I8-1) / .381830050505119E0,
+ > .279705391489277E0,
+ > .129484966168870E0/
+C N=8
+ DATA (RZKSI(I),I=I8,I9-1) / .183434642495650E0,
+ > .525532409916329E0,
+ > .796666477413627E0,
+ > .960289856497536E0/
+ DATA (RWKSI(I),I=I8,I9-1) / .362683783378362E0,
+ > .313706645877887E0,
+ > .222381034453374E0,
+ > .101228536290376E0/
+C N=9
+ DATA (RZKSI(I),I=I9,I10-1) / .324253423403809E0,
+ > .613371432700590E0,
+ > .836031107326636E0,
+ > .968160239507626E0/
+ DATA (RWKSI(I),I=I9,I10-1) / .312347077040003E0,
+ > .260610696402935E0,
+ > .180648160694857E0,
+ > .081274388361574E0/
+C N=10
+ DATA (RZKSI(I),I=I10,I11-1) / .148874338981631E0,
+ > .433395394129247E0,
+ > .679409568299024E0,
+ > .865063366688985E0,
+ > .973906528517172E0/
+ DATA (RWKSI(I),I=I10,I11-1) / .295524224714753E0,
+ > .269266719309996E0,
+ > .219086362515982E0,
+ > .149451349150581E0,
+ > .066671344308688E0/
+C N=11
+ DATA (RZKSI(I),I=I11,I12-1) /
+ > .978228658146057E0,.887062599768095E0,
+ > .730152005574049E0,.519096129206812E0,
+ > .269543155952345E0/
+ DATA (RWKSI(I),I=I11,I12-1) /
+ > .055668567116174E0,.125580369464905E0,
+ > .186290210927734E0,.233193764591990E0,
+ > .262804544510247E0/
+C N=12
+ DATA (RZKSI(I),I=I12,I13-1) /
+ > .125233408511469E0,.367831498998180E0,
+ > .587317954286617E0,.769902674194305E0,
+ > .904117256370475E0,.981560634246719E0/
+ DATA (RWKSI(I),I=I12,I13-1) /
+ > .249147045813403E0,.233492536538355E0,
+ > .203167426723066E0,.160078328543346E0,
+ > .106939325995318E0,.047175336386512E0/
+C N=13
+ DATA (RZKSI(I),I=I13,I14-1) /
+ > .984183054718588E0,.917598399222978E0,
+ > .801578090733310E0,.642349339440340E0,
+ > .448492751036447E0,.230458315955135E0/
+ DATA (RWKSI(I),I=I13,I14-1) /
+ > .040484004765316E0,.092121499837728E0,
+ > .138873510219787E0,.178145980761946E0,
+ > .207816047536888E0,.226283180262897E0/
+C N=14
+ DATA (RZKSI(I),I=I14,I15-1) /
+ > .986283808696812E0,.928434883663573E0,
+ > .827201315069765E0,.687292904811685E0,
+ > .515248636358154E0,.319112368927890E0,
+ > .108054948707344E0/
+ DATA (RWKSI(I),I=I14,I15-1) /
+ > .035119460331752E0,.080158087159760E0,
+ > .121518570687903E0,.157203167158193E0,
+ > .185538397477938E0,.205198463721296E0,
+ > .215263853463158E0/
+C N=15
+ DATA (RZKSI(I),I=I15,I16-1) /
+ > .987992518020485E0,.937273392400706E0,
+ > .848206583410427E0,.724417731360170E0,
+ > .570972172608539E0,.394151347077563E0,
+ > .201194093997434E0/
+ DATA (RWKSI(I),I=I15,I16-1) /
+ > .030753241996117E0,.070366047488108E0,
+ > .107159220467172E0,.139570677926154E0,
+ > .166269205816994E0,.186161000015562E0,
+ > .198431485327111E0/
+C N=16
+ DATA (RZKSI(I),I=I16,I17-1) /
+ > .095012509837637E0,.281603550779259E0,
+ > .458016777657227E0,.617876244402644E0,
+ > .755404408355003E0,.865631202387832E0,
+ > .944575023073233E0,.989400934991650E0/
+ DATA (RWKSI(I),I=I16,I17-1) /
+ > .189450610455068E0,.182603415044924E0,
+ > .169156519395003E0,.149595988816577E0,
+ > .124628971255534E0,.095158511682493E0,
+ > .062253523938648E0,.027152459411754E0/
+C N=17
+ DATA (RZKSI(I),I=I17,I18-1) /
+ > .990575475314417E0,.950675521768768E0,
+ > .880239153726986E0,.781514003896801E0,
+ > .657671159216691E0,.512690537086477E0,
+ > .351231763453876E0,.178484181495848E0/
+ DATA (RWKSI(I),I=I17,I18-1) /
+ > .024148302868548E0,.055459529373987E0,
+ > .085036148317179E0,.111883847193404E0,
+ > .135136368468525E0,.154045761076810E0,
+ > .168004102156450E0,.176562705366993E0/
+C N=18
+ DATA (RZKSI(I),I=I18,I19-1) /
+ > .991565168420930E0,.955823949571397E0,
+ > .892602466497556E0,.803704958972523E0,
+ > .691687043060353E0,.559770831073947E0,
+ > .411751161462843E0,.251886225691505E0,
+ > .084775013041735E0/
+ DATA (RWKSI(I),I=I18,I19-1) /
+ > .021616013526483E0,.049714548894970E0,
+ > .076425730254889E0,.100942044106287E0,
+ > .122555206711478E0,.140642914670651E0,
+ > .154684675126265E0,.164276483745833E0,
+ > .169142382963143E0/
+C N=19
+ DATA (RZKSI(I),I=I19,I20-1) /
+ > .992406843843584E0,.960208152134830E0,
+ > .903155903614818E0,.822714656537143E0,
+ > .720966177335229E0,.600545304661681E0,
+ > .464570741375961E0,.316564099963630E0,
+ > .160358645640225E0/
+ DATA (RWKSI(I),I=I19,I20-1) /
+ > .019461788229726E0,.044814226765700E0,
+ > .069044542737641E0,.091490021622450E0,
+ > .111566645547334E0,.128753962539336E0,
+ > .142606702173607E0,.152766042065860E0,
+ > .158968843393954E0/
+C N=20
+ DATA (RZKSI(I),I=I20,I24-1) /
+ > .076526521133497E0,.227785851141645E0,
+ > .373706088715420E0,.510867001950827E0,
+ > .636053680726515E0,.746331906460151E0,
+ > .839116971822219E0,.912234428251326E0,
+ > .963971927277914E0,.993128599185095E0/
+ DATA (RWKSI(I),I=I20,I24-1) /
+ > .152753387130726E0,.149172986472604E0,
+ > .142096109318382E0,.131688638449177E0,
+ > .118194531961518E0,.101930119817240E0,
+ > .083276741576705E0,.062672048334109E0,
+ > .040601429800387E0,.017614007139152E0/
+C N=24
+ DATA (RZKSI(I),I=I24,I28-1) /
+ > .995187219997021E0,.974728555971309E0,
+ > .938274552002733E0,.886415527004401E0,
+ > .820001985973903E0,.740124191578554E0,
+ > .648093651936975E0,.545421471388839E0,
+ > .433793507626045E0,.315042679696163E0,
+ > .191118867473616E0,.064056892862605E0/
+ DATA (RWKSI(I),I=I24,I28-1) /
+ > .012341229799987E0,.028531388628934E0,
+ > .044277438817420E0,.059298584915437E0,
+ > .073346481411080E0,.086190161531953E0,
+ > .097618652104114E0,.107444270115966E0,
+ > .115505668053726E0,.121670472927803E0,
+ > .125837456346828E0,.127938195346752E0/
+C N=28
+ DATA (RZKSI(I),I=I28,I32-1) /
+ > .996442497573954E0,.981303165370873E0,
+ > .954259280628938E0,.915633026392132E0,
+ > .865892522574395E0,.805641370917179E0,
+ > .735610878013632E0,.656651094038865E0,
+ > .569720471811402E0,.475874224955118E0,
+ > .376251516089079E0,.272061627635178E0,
+ > .164569282133380E0,.055079289884034E0/
+ DATA (RWKSI(I),I=I28,I32-1) /
+ > .009124282593094E0,.021132112592771E0,
+ > .032901427782304E0,.044272934759004E0,
+ > .055107345675717E0,.065272923966999E0,
+ > .074646214234569E0,.083113417228901E0,
+ > .090571744393033E0,.096930657997930E0,
+ > .102112967578061E0,.106055765922846E0,
+ > .108711192258294E0,.110047013016475E0/
+C N=32
+ DATA (RZKSI(I),I=I32,I64-1) /
+ > .048307665687738E0,.144471961582796E0,
+ > .239287362252137E0,.331868602282128E0,
+ > .421351276130635E0,.506899908932229E0,
+ > .587715757240762E0,.663044266930215E0,
+ > .732182118740290E0,.794483795967942E0,
+ > .849367613732570E0,.896321155766052E0,
+ > .934906075937740E0,.964762255587506E0,
+ > .985611511545268E0,.997263861849482E0/
+ DATA (RWKSI(I),I=I32,I64-1) /
+ > .096540088514728E0,.095638720079275E0,
+ > .093844399080805E0,.091173878695764E0,
+ > .087652093004404E0,.083311924226947E0,
+ > .078193895787070E0,.072345794108849E0,
+ > .065822222776362E0,.058684093478536E0,
+ > .050998059262376E0,.042835898022227E0,
+ > .034273862913021E0,.025392065309262E0,
+ > .016274394730906E0,.007018610009470E0/
+C N=64
+ DATA (RZKSI(I),I=I64,IFN) /
+ > .024350292663424E0,.072993121787799E0,.121462819296121E0,
+ > .169644420423993E0,.217423643740007E0,.264687162208767E0,
+ > .311322871990211E0,.357220158337668E0,.402270157963992E0,
+ > .446366017253464E0,.489403145707053E0,.531279464019894E0,
+ > .571895646202634E0,.611155355172393E0,.648965471254657E0,
+ > .685236313054233E0,.719881850171611E0,.752819907260532E0,
+ > .783972358943341E0,.813265315122798E0,.840629296252580E0,
+ > .865999398154093E0,.889315445995114E0,.910522137078503E0,
+ > .929569172131940E0,.946411374858403E0,.961008799652054E0,
+ > .973326827789911E0,.983336253884626E0,.991013371476744E0,
+ > .996340116771955E0,.999305041735772E0/
+ DATA (RWKSI(I),I=I64,IFN) /
+ > .048690957009140E0,.048575467441503E0,.048344762234803E0,
+ > .047999388596458E0,.047540165714830E0,.046968182816210E0,
+ > .046284796581314E0,.045491627927418E0,.044590558163757E0,
+ > .043583724529323E0,.042473515123654E0,.041262563242624E0,
+ > .039953741132720E0,.038550153178616E0,.037055128540240E0,
+ > .035472213256882E0,.033805161837142E0,.032057928354852E0,
+ > .030234657072402E0,.028339672614259E0,.026377469715055E0,
+ > .024352702568711E0,.022270173808383E0,.020134823153530E0,
+ > .017951715775697E0,.015726030476025E0,.013463047896719E0,
+ > .011168139460131E0,.008846759826364E0,.006504457968978E0,
+ > .004147033260562E0,.001783280721696E0/
+*
+ IDEP=0
+ IFIN=0
+ FW=HALF*(XSUP-XINF)
+ FC=HALF*(XSUP+XINF)
+ IF( NGPT.EQ. 1 ) THEN
+ ZGKSI(1)=FC
+ WGKSI(1)=TWO*FW
+ RETURN
+ ELSE IF( NGPT.EQ. 2 ) THEN
+ IDEP=I2
+ IFIN=I3-1
+ ELSE IF( NGPT.EQ. 3 ) THEN
+ IDEP=I3
+ IFIN=I4-1
+ ELSE IF( NGPT.EQ. 4 ) THEN
+ IDEP=I4
+ IFIN=I5-1
+ ELSE IF( NGPT.EQ. 5 ) THEN
+ IDEP=I5
+ IFIN=I6-1
+ ELSE IF( NGPT.EQ. 6 ) THEN
+ IDEP=I6
+ IFIN=I7-1
+ ELSE IF( NGPT.EQ. 7 ) THEN
+ IDEP=I7
+ IFIN=I8-1
+ ELSE IF( NGPT.EQ. 8 ) THEN
+ IDEP=I8
+ IFIN=I9-1
+ ELSE IF( NGPT.EQ. 9 ) THEN
+ IDEP=I9
+ IFIN=I10-1
+ ELSE IF( NGPT.EQ.10 ) THEN
+ IDEP=I10
+ IFIN=I11-1
+ ELSE IF( NGPT.EQ.11 ) THEN
+ IDEP=I11
+ IFIN=I12-1
+ ELSE IF( NGPT.EQ.12 ) THEN
+ IDEP=I12
+ IFIN=I13-1
+ ELSE IF( NGPT.EQ.13 ) THEN
+ IDEP=I13
+ IFIN=I14-1
+ ELSE IF(NGPT.EQ.14 ) THEN
+ IDEP=I14
+ IFIN=I15-1
+ ELSE IF(NGPT.EQ.15 ) THEN
+ IDEP=I15
+ IFIN=I16-1
+ ELSE IF(NGPT.EQ.16 ) THEN
+ IDEP=I16
+ IFIN=I17-1
+ ELSE IF(NGPT.EQ.17 ) THEN
+ IDEP=I17
+ IFIN=I18-1
+ ELSE IF(NGPT.EQ.18 ) THEN
+ IDEP=I18
+ IFIN=I19-1
+ ELSE IF(NGPT.EQ.19 ) THEN
+ IDEP=I19
+ IFIN=I20-1
+ ELSE IF( NGPT.EQ.20 ) THEN
+ IDEP=I20
+ IFIN=I24-1
+ ELSE IF( NGPT.EQ.24 ) THEN
+ IDEP=I24
+ IFIN=I28-1
+ ELSE IF( NGPT.EQ.28 ) THEN
+ IDEP=I28
+ IFIN=I32-1
+ ELSE IF( NGPT.EQ.32 ) THEN
+ IDEP=I32
+ IFIN=I64-1
+ ELSE IF( NGPT.EQ.64 ) THEN
+ IDEP=I64
+ IFIN=IFN
+ ELSE
+ WRITE(CERROR,'(I4)') NGPT
+ CALL XABORT('ALGPT:INVALID NBR.OF GAUSS PTS.='//CERROR//' ***
+ > 1 TO 20,24,28,32,64 PTS. ARE PERMITTED')
+ ENDIF
+C------
+C INITIALIZE ZGKSI AND WGKSI STARTING FROM BOTTOM WITH NEGATIVE
+C AND FINISHING WITH POSITIVE
+C------
+ IUP=1
+ IDN=NGPT
+ DO 100 I=IFIN,IDEP,-1
+ ZGKSI(IUP)=-FW*RZKSI(I)+FC
+ WGKSI(IUP)=FW*RWKSI(I)
+ ZGKSI(IDN)=FW*RZKSI(I)+FC
+ WGKSI(IDN)=FW*RWKSI(I)
+ IUP=IUP+1
+ IDN=IDN-1
+ 100 CONTINUE
+C------
+C FOR ODD NGPT, CENTRAL POINT IS 0.0 AND WEIGHT IS (2.0-TOTAL WEIGHT)
+C------
+ IF(IUP.EQ.IDN) THEN
+ ZGKSI(IUP)=FC
+ DO 110 I=1,IUP-1
+ FW=FW-WGKSI(I)
+ 110 CONTINUE
+ WGKSI(IUP)=TWO*FW
+ ENDIF
+ RETURN
+ END