path: root/Dragon/src/SNQU05.f

*DECK SNQU05
      SUBROUTINE SNQU05(NLF,TPQ,UPQ,VPQ,WPQ)
*
*-----------------------------------------------------------------------
*
*Purpose:
* To define symmetric Legendre-Chebyshev quadrature angles.
*
*Copyright:
* Copyright (C) 2005 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): 
* M. Hampartzounian 
*
*Update(s): 
* Adapted to SNT: module of N.Martin.
*
*Parameters: input
* NLF  order of the SN approximation (even number). 
*
*Parameters: output
* TPQ     base points in $\\xi$ of the 2D/3D SN quadrature.
* UPQ     base points in $\\mu$ of the 2D/3D SN quadrature.
* VPQ     base points in $\\eta$ of the 2D/3D SN quadrature.
* WPQ     weights of the 2D/3D SN quadrature.
*
*Reference: 
* G. Longoni, A. Haghighat, 
* Development of New Quadrature Sets with the 
* Ordinate Splitting Technique,
* M&C 2001, Salt Lake City, September 2001.  
*
*-----------------------------------------------------------------------
*
*----
*  Subroutine arguments
*----
      INTEGER NLF
      REAL TPQ(NLF*(NLF/2+1)/4),UPQ(NLF*(NLF/2+1)/4),
     1     VPQ(NLF*(NLF/2+1)/4),WPQ(NLF*(NLF/2+1)/4)
*----
*  Local parameters
*----
      INTEGER          IOUT
      CHARACTER        NAMSBR*6
      PARAMETER       (IOUT=6,NAMSBR='SNQU05',MAXNLF=44)
      DOUBLE PRECISION DZERO,DONE,DTWO,PID2
      PARAMETER       (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
*----
*  Functions
*----
      DOUBLE PRECISION XDRCST,PI
*----
*  Local variables
*----
      INTEGER          NBPT,ILEVEL,JLEVEL,NLEVEL,KLEVEL,IPOINT,IPT
      DOUBLE PRECISION ALPHA,XI,ROTXI,COSA,SINA
*----
*  Data
*----
      INTEGER          IS,IA
      DOUBLE PRECISION SYST(2,253)
      SAVE             SYST
*----
*  Set quadrature points and weights
*----
      IF(NLF.GT.MAXNLF) CALL XABORT('SNQU05: MAXNLF OVERFLOW.')
*----
*  Order =   2
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=  1,  1)/ 
     > 0.577350269189626D0,1.0D0/
*----
*  Order =   4
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=  2,  3)/
     > 0.33998104358486D0,0.32607257743127D0,0.86113631159405D0,
     > 0.34785484513745D0/
*----
*  Order =   6
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=  4,  6)/
     > 0.23861918608320D0,0.15597131152423D0,0.66120938646626D0,
     > 0.18038078652407D0,0.93246951420315D0,0.17132449237917D0/
*----
*  Order =   8
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=  7, 10)/
     > 0.18343464249565D0,0.09067094584459D0,0.52553240991633D0,
     > 0.10456888195930D0,0.79666647741362D0,0.11119051722669D0,
     > 0.96028985649754D0,0.10122853629037D0/
*----
*  Order =  10
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 11, 15)/
     > 0.14887433898163D0,0.05910484494295D0,0.43339539412925D0,
     > 0.06731667982750D0,0.67940956829902D0,0.07302878750533D0,
     > 0.86506336668898D0,0.07472567457529D0,0.97390652851718D0,
     > 0.06667134430868D0/
*----
*  Order =  12
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 16, 21)/
     > 0.12523340851147D0,0.04152450763557D0,0.36783149899818D0,
     > 0.04669850730767D0,0.58731795428663D0,0.05079185668077D0,
     > 0.76990267419428D0,0.05335944284779D0,0.90411725637050D0,
     > 0.05346966299765D0,0.98156063424671D0,0.04717533638653D0/
*----
*  Order =  14
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 22, 28)/
     > 0.10805494870734D0,0.03075197906617D0,0.31911236892789D0,
     > 0.03419974395355D0,0.51524863635816D0,0.03710767949559D0,
     > 0.68729290481172D0,0.03930079178955D0,0.82720131506965D0,
     > 0.04050619022934D0,0.92843488366371D0,0.04007904357982D0,
     > 0.98628380869675D0,0.03511946033190D0/
*----
*  Order =   16
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 29, 36)/
     > 0.09501250983764D0,0.02368132630688D0,0.28160355077926D0,
     > 0.02608620214927D0,0.45801677765726D0,0.02819275323250D0,
     > 0.61787624440245D0,0.02991919776332D0,0.75540440835553D0,
     > 0.03115724281390D0,0.86563120238710D0,0.03171950389411D0,
     > 0.94457502307371D0,0.03112676196937D0,0.98940093499154D0,
     > 0.02715245941196D0/
*----
*  Order =   18
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 37, 45)/
     > 0.08477501304174D0,0.01879359810702D0,0.25188622569150D0,
     > 0.02053456046823D0,0.41175116146291D0,0.02209781073233D0,
     > 0.55977083107372D0,0.02344048577842D0,0.69168704306110D0,
     > 0.02451104134236D0,0.80370495897046D0,0.02523551102634D0,
     > 0.89260246650148D0,0.02547524341783D0,0.95582394956707D0,
     > 0.02485727444991D0,0.99156516842283D0,0.02161601352180D0/
*----
*  Order =  20
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 46, 55)/
     > 0.07652652113350D0,0.01527533871307D0,0.22778585114165D0,
     > 0.01657477627473D0,0.37370608871523D0,0.01776201366480D0,
     > 0.51086700195199D0,0.01881266263554D0,0.63605368072293D0,
     > 0.01969908866072D0,0.74633190646642D0,0.02038602396117D0,
     > 0.83911697181714D0,0.02081918540086D0,0.91223442824898D0,
     > 0.02089068276635D0,0.96397192728679D0,0.02030071490687D0,
     > 0.99312859917998D0,0.01761400714924D0/
*----
*  Order =  22
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 56, 66)/
     > 0.06973927331972D0,0.01265926116869D0,0.20786042668821D0,
     > 0.01365414983460D0,0.34193582089230D0,0.01457483386523D0,
     > 0.46935583798525D0,0.01540654710137D0,0.58764040351215D0,
     > 0.01613318515391D0,0.69448726317471D0,0.01673569074235D0,
     > 0.78781680600113D0,0.01718832124013D0,0.86581257768321D0,
     > 0.01744911712470D0,0.92695677224150D0,0.01743111170530D0,
     > 0.97006049778158D0,0.01688745083050D0,0.99429458550512D0,
     > 0.01462799524418D0/
*----
*  Order =  24
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 67, 78)/
     > 0.06405689286261D0,0.01066151627890D0,0.19111886747361D0,
     > 0.01143976875880D0,0.31504267969640D0,0.01216704729278D0,
     > 0.43379350762497D0,0.01283396311708D0,0.54542147138999D0,
     > 0.01343053376483D0,0.64809365193896D0,0.01394552172544D0,
     > 0.74012419158451D0,0.01436502694016D0,0.82000198588723D0,
     > 0.01466929621748D0,0.88641552731297D0,0.01482464636908D0,
     > 0.93827455142585D0,0.01475914616380D0,0.97472855656607D0,
     > 0.01426569422080D0,0.99518721974900D0,0.01234123034267D0/
*----
*  Order =  26
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 79, 91)/
     > 0.05923009342931D0,0.00910164732917D0,0.17685882035693D0,
     > 0.00972170362377D0,0.29200483948558D0,0.01030561968603D0,
     > 0.40305175512448D0,0.01084718405284D0,0.50844071482581D0,
     > 0.01133990678847D0,0.60669229300123D0,0.01177672504332D0,
     > 0.69642726048349D0,0.01214941348255D0,0.77638594865137D0,
     > 0.01244735829575D0,0.84544594311257D0,0.01265480931773D0,
     > 0.90263786157297D0,0.01274395622539D0,0.94715906696364D0,
     > 0.01265412770351D0,0.97838544586603D0,0.01220892547636D0,
     > 0.99588570114160D0,0.01055137268743D0/
*----
*  Order =  28
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA= 92,105)/
     > 0.05507928988403D0,0.00786050092975D0,0.16456928213333D0,
     > 0.00836239940448D0,0.27206162763654D0,0.00883798049356D0,
     > 0.37625151607709D0,0.00928299705271D0,0.47587422501394D0,
     > 0.00969306579864D0,0.56972047161791D0,0.01006352716248D0,
     > 0.65665109450180D0,0.01038917711347D0,0.73561087718791D0,
     > 0.01066374513244D0,0.80564137195850D0,0.01087881964355D0,
     > 0.86589252190936D0,0.01102147166953D0,0.91563302581684D0,
     > 0.01106822904179D0,0.95425928269356D0,0.01096714930000D0,
     > 0.98130316299531D0,0.01056605251579D0,0.99644249859202D0,
     > 0.00912428140132D0/
*----
*  Order =  30
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=106,120)/
     > 0.05147184255532D0,0.00685684352624D0,0.15386991360868D0,
     > 0.00726874212489D0,0.25463692616461D0,0.00766103235284D0,
     > 0.35270472554884D0,0.00803072809779D0,0.44703376949342D0,
     > 0.00837477474924D0,0.53662414821306D0,0.00868997871608D0,
     > 0.62052618288637D0,0.00897287728847D0,0.69785049482498D0,
     > 0.00921949657093D0,0.76777743319354D0,0.00942489107957D0,
     > 0.82956575621336D0,0.00958219000643D0,0.88256055526296D0,
     > 0.00968054393272D0,0.92620000651698D0,0.00969978543893D0,
     > 0.96002192312100D0,0.00959491170268D0,0.98366807161209D0,
     > 0.00923324554691D0,0.99689350414176D0,0.00796814638786D0/
*----
*  Order =  32
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=121,136)/
     > 0.04830766568774D0,0.00603375553217D0,0.14447196158186D0,
     > 0.00637591467195D0,0.23928736227279D0,0.00670317136284D0,
     > 0.33186860213190D0,0.00701337528509D0,0.42135127674309D0,
     > 0.00730434107922D0,0.50689990719947D0,0.00757381130805D0,
     > 0.58771576112490D0,0.00781938956129D0,0.66304425908752D0,
     > 0.00803842135657D0,0.73218213532892D0,0.00822777972086D0,
     > 0.79448375640681D0,0.00838343484386D0,0.84936770896391D0,
     > 0.00849970145691D0,0.89632095843854D0,0.00856713191057D0,
     > 0.93490639505449D0,0.00856850293792D0,0.96476187795073D0,
     > 0.00846401442985D0,0.98561180730897D0,0.00813711928812D0,
     > 0.99726375488382D0,0.00701885042221D0/
*----
*  Order =  34
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=137,153)/
     > 0.04550982195310D0,0.00535039649002D0,0.13615235725880D0,
     > 0.00563769027317D0,0.22566669163722D0,0.00591345985566D0,
     > 0.31331108104631D0,0.00617612426888D0,0.39835927975331D0,
     > 0.00642408458066D0,0.48010653651464D0,0.00665570376058D0,
     > 0.55787552837148D0,0.00686927018743D0,0.63102165671823D0,
     > 0.00706293794571D0,0.69893926435390D0,0.00723461457806D0,
     > 0.76106458437911D0,0.00738175933892D0,0.81688476588987D0,
     > 0.00750110078301D0,0.86593368730159D0,0.00758751425230D0,
     > 0.90781120740225D0,0.00763351850909D0,0.94216032905008D0,
     > 0.00762273892258D0,0.96871043413727D0,0.00752113249273D0,
     > 0.98722623909002D0,0.00722576516252D0,0.99757230194591D0,
     > 0.00622779370687D0/
*----
*  Order =  36
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=154,171)/
     > 0.04301819847371D0,0.00477684864836D0,0.12873610380802D0,
     > 0.00502039327879D0,0.21350089236042D0,0.00525488868613D0,
     > 0.29668499489157D0,0.00547915111524D0,0.37767254966703D0,
     > 0.00569198776656D0,0.45586393486895D0,0.00589218549996D0,
     > 0.53068031264856D0,0.00607849018656D0,0.60156759861950D0,
     > 0.00624957522881D0,0.66800134713266D0,0.00640398037611D0,
     > 0.72948900231201D0,0.00654000932504D0,0.78557638994910D0,
     > 0.00665561122424D0,0.83584737503067D0,0.00674782531543D0,
     > 0.87992821829261D0,0.00681274160010D0,0.91750240344171D0,
     > 0.00684271263501D0,0.94826416091557D0,0.00682558556827D0,
     > 0.97203921707055D0,0.00672476115101D0,0.98857686223253D0,
     > 0.00646291222847D0,0.99783406138946D0,0.00555667184560D0/
*----
*  Order =  38
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=172,190)/
     > 0.04078514790459D0,0.00429079101476D0,0.12208402533452D0,
     > 0.00449902743170D0,0.20257045398628D0,0.00470006077885D0,
     > 0.28170880888511D0,0.00489299029385D0,0.35897244495443D0,
     > 0.00507691088444D0,0.43384715631935D0,0.00525090669126D0,
     > 0.50583473801940D0,0.00541403881989D0,0.57445603748646D0,
     > 0.00556532827760D0,0.63925418826503D0,0.00570372124187D0,
     > 0.69979962680485D0,0.00582803437120D0,0.75568295130273D0,
     > 0.00593692197056D0,0.80655201425815D0,0.00602847883995D0,
     > 0.85201700089545D0,0.00610101854081D0,0.89189086124464D0,
     > 0.00614733110727D0,0.92568500732863D0,0.00617630505125D0,
     > 0.95353898729968D0,0.00613319286783D0,0.97477439677976D0,
     > 0.00607355513150D0,0.98978922946055D0,0.00577832504734D0,
     > 0.99803301416225D0,0.00504805879326D0/
*----
*  Order =  40
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=191,210)/
     > 0.03877241750606D0,0.00387529739892D0,0.11608407067269D0,
     > 0.00405472727181D0,0.19269758075384D0,0.00422835343883D0,
     > 0.26815218442225D0,0.00439548053431D0,0.34199409501760D0,
     > 0.00455541138509D0,0.41377918521708D0,0.00470744324502D0,
     > 0.48307585742886D0,0.00485086018626D0,0.54946703621944D0,
     > 0.00498492376041D0,0.61255384817576D0,0.00510886273005D0,
     > 0.67195746920806D0,0.00522174107642D0,0.72731518674779D0,
     > 0.00532302904828D0,0.77831437316505D0,0.00540964292729D0,
     > 0.82459057917465D0,0.00548656553689D0,0.86600919232031D0,
     > 0.00553340517150D0,0.90199639822285D0,0.00558997462627D0,
     > 0.93299517738927D0,0.00556748623087D0,0.95765721118329D0,
     > 0.00561397264950D0,0.97754603471325D0,0.00541857778427D0,
     > 0.99051106402526D0,0.00539644941610D0,0.99831206210441D0,
     > 0.00434333104876D0/
*----
*  Order =  42
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=211,231)/
     > 0.03694894316534D0,0.00351734448725D0,0.11064502720743D0,
     > 0.00367304067267D0,0.18373680686444D0,0.00382400922292D0,
     > 0.25582507379698D0,0.00396970636014D0,0.32651616975104D0,
     > 0.00410958765725D0,0.39542362251945D0,0.00424310738074D0,
     > 0.46217274386666D0,0.00436970393775D0,0.52639343800021D0,
     > 0.00448884199956D0,0.58774968956786D0,0.00459981369037D0,
     > 0.64587456004057D0,0.00470237199319D0,0.70050533824065D0,
     > 0.00479452034217D0,0.75127821879747D0,0.00488025174105D0,
     > 0.79791981602651D0,0.00494337495192D0,0.84045950845072D0,
     > 0.00502154054607D0,0.87755806345153D0,0.00503689455872D0,
     > 0.91197509273697D0,0.00517770715670D0,0.93729218041353D0,
     > 0.00526592418727D0,0.96414624148124D0,0.00570793259377D0,
     > 0.97695704453476D0,0.00690539560442D0,0.99322113425436D0,
     > 0.00541384107796D0,0.99769796354612D0,0.00754180253610D0/
*----
*  Order =  44
*----
      DATA             ((SYST(IS,IA),IS=1,2),IA=232,253)/
     > 0.03528923696410D0,0.00320677989952D0,0.10569190169021D0,
     > 0.00334274692733D0,0.17556801710226D0,0.00347482459015D0,
     > 0.24456941646933D0,0.00360258272163D0,0.31235278934078D0,
     > 0.00372559021852D0,0.37857772181427D0,0.00384342407853D0,
     > 0.44292618820931D0,0.00395560422690D0,0.50503688871692D0,
     > 0.00406188071758D0,0.56471466385375D0,0.00416116332728D0,
     > 0.62139067754562D0,0.00425481282618D0,0.67534196136471D0,
     > 0.00433842941591D0,0.72528399895512D0,0.00441848009213D0,
     > 0.77261708065981D0,0.00448920279793D0,0.81468343050057D0,
     > 0.00453553695561D0,0.85451710591231D0,0.00462393394081D0,
     > 0.88796251449581D0,0.00457908321421D0,0.91913386326312D0,
     > 0.00479985160516D0,0.94407409289036D0,0.00460697658507D0,
     > 0.96484471190689D0,0.00431337311996D0,0.98174671963958D0,
     > 0.00476874072667D0,0.99175403861920D0,0.00484009055849D0,
     > 0.99874459587735D0,0.00339611207931D0/
*----
*  Start processing
*----        
      PI=XDRCST('Pi',' ')
      PID2=PI/DTWO
      NBPT=(NLF*(NLF+2))/8
      ILEVEL=1
      JLEVEL=0
      NLEVEL=NLF/2
      KLEVEL=(NLEVEL*(NLEVEL-1))/2
      IPOINT=0
      DO IPT=1,NBPT
        IPOINT=IPOINT+1
        JLEVEL = JLEVEL + 1
        ALPHA =PI*(DBLE(NLEVEL-ILEVEL+1-JLEVEL)+DONE/DTWO)
     >         /(DTWO*DBLE(NLEVEL-ILEVEL+1))
        XI=SYST(1,KLEVEL+ILEVEL)
        ROTXI=SQRT(DONE-XI**2)
        COSA=COS(ALPHA)
        SINA=SIN(ALPHA)
        TPQ(IPOINT)=REAL(COSA*ROTXI)
        UPQ(IPOINT)=REAL(SINA*ROTXI)
        VPQ(IPOINT)=REAL(XI)
        WPQ(IPOINT)=REAL(PID2*SYST(2,KLEVEL+ILEVEL))
        IF(JLEVEL .EQ. NLEVEL-ILEVEL+1) THEN 
          ILEVEL = ILEVEL + 1
           JLEVEL = 0
        ENDIF
      ENDDO
      RETURN
      END