separatrice

PURPOSE ^

Fortran source

SYNOPSIS ^

Fortran source

DESCRIPTION ^

Fortran source

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 C
0002 #include "fintrf.h"
0003 C
0004       subroutine mexFunction(nlhs,plhs,nrhs,prhs)
0005 C
0006 C Compilation : !mex -f mexopts_f77.sh separatrice.f (dans matlab7)
0007 C--------------------------------------------------------------
0008 
0009       IMPLICIT REAL*8 (A-H,O-Z)
0010       IMPLICIT integer (I-N)
0011 C DECLARATION DES VARIABLES
0012 c pointeurs mex   
0013       mwPointer plhs(*), prhs(*)
0014 
0015       integer nlhs, nrhs
0016 
0017       mwPointer mxGetPr, mxCreateDoubleMatrix
0018 
0019       integer mxgetN, mxgetM
0020 c 
0021 c Extrait du main de helios.f concernant le trace de la derniere surface
0022 c magnetique (il s'agit d'une implementation triviale des formules fournies)
0023 c
0024 c pointeurs et tailles   
0025    
0026       mwPointer geom_pr,mode_pr,x_pr,y_pr
0027       integer ifail,iverbose
0028       real*8 mode,geom(8),x(500),y(500)
0029 
0030       mwSize mat1,mat8,mat500
0031       parameter (mat1=1) 
0032       parameter (mat8=8)
0033       parameter (mat500=500)
0034       
0035       character(len=200) msg
0036 
0037       data pi/3.141592653589793/
0038       data two/2./,one/1./,rho/1./ 
0039 c      data rkappax1/1.687/,ptdeltax1/0.466/,psip1dg/0./,psim1dg/0./,
0040 c     & rkappax2/2.001/,ptdeltax2/0.568/,psip2dg/22.46/,psim2dg/67.92/
0041 **********************************************************************
0042 
0043 c verbose mode 
0044 
0045      iverbose = 0
0046 
0047 C---------------------------------------------------------------------
0048 C GESTION DES ENTREES     
0049 C entree
0050        
0051       if(iverbose .eq. 1)then
0052       write(msg,*) "debut separatrice\n" 
0053       call mexPrintf(msg)
0054       endif  
0055       geom_pr = mxGetPr(prhs(1))
0056       call mxCopyPtrToReal8(geom_pr,geom,mat8)
0057  
0058       rkappax1  = geom(1)
0059       ptdeltax1 = geom(2)
0060       psip1dg   = geom(3)
0061       psim1dg   = geom(4)
0062       rkappax2  = geom(5)
0063       ptdeltax2 = geom(6)
0064       psip2dg   = geom(7)
0065       psim2dg   = geom(8)
0066 C mode de calcul
0067       mode_pr = mxGetPr(prhs(2))
0068       call mxCopyPtrToReal8(mode_pr,mode,mat1)
0069 
0070 c Initializing geometrical parameters
0071       psimoins1=psim1dg*pi/180.
0072       psiplus1=psip1dg*pi/180.
0073       psimoins2=psim2dg*pi/180.
0074       psiplus2=psip2dg*pi/180.
0075       tmoins1=tan(psimoins1)*(1.-ptdeltax1)/rkappax1
0076       tmoins2=tan(psimoins2)*(1.-ptdeltax2)/rkappax2
0077       if(iverbose .eq. 1)then
0078       write(msg,*) 'tmoins1 = ',tmoins1,' tmoins2 =',tmoins2,"\n" 
0079       call mexPrintf(msg)
0080       endif
0081       tplus1=tan(psiplus1)*(1.+ptdeltax1)/rkappax1
0082       tplus2=tan(psiplus2)*(1.+ptdeltax2)/rkappax2
0083       if(iverbose .eq. 1)then
0084       write(msg,*) "tplus1 = ",tplus1," tplus2 =",tplus2,"\n"
0085       call mexPrintf(msg)
0086       endif
0087 c      
0088 c Activer manuellement la bonne partie du programme suivant le cas (voir texte
0089 c fourni) en remplacant if(one.gt.two) par if(two.gt.one)
0090 c Un algorithme automatique peut aussi 
0091 c
0092 c plot of the plasma last closed magnetic surface
0093 c
0094 c case where the upper part is made of 2 portions of ellipse and the
0095 c bottom part is made of 2 portions of ellipse
0096 c
0097       if(mode .eq. 1)then
0098         tetaxp2   = 0.5*pi-asin(tplus2/(1.-tplus2))
0099         tetap2in  = 0.
0100         tetap2fi  = tetaxp2
0101         nbtetap2  = 125
0102         tetap2st  = (tetap2fi-tetap2in)/(nbtetap2-1)
0103         tetaxm2   = 0.5*pi+asin(tmoins2/(1.-tmoins2))
0104         tetam2in  = tetaxm2
0105         tetam2fi  = pi
0106         nbtetam2  = 125
0107         tetam2st  = (tetam2fi-tetam2in)/(nbtetam2-1)
0108         tetaxm1   = 0.5*pi+asin(tmoins1/(1.-tmoins1))
0109         tetam1in  = pi
0110         tetam1fi  = tetaxm1
0111         nbtetam1  = 125
0112         tetam1st  = (tetam1fi-tetam1in)/(nbtetam1-1)
0113         tetaxp1   = 0.5*pi-asin(tplus1/(1.-tplus1))
0114         tetap1in  = -tetaxp1
0115         tetap1fi  = 0.
0116         nbtetap1  = 125
0117         tetap1st  = (tetap1fi-tetap1in)/(nbtetap1-1)
0118         ii        = 1
0119 c bottom plus
0120         do noteta=1,nbtetap2
0121           teta     = tetap2in+(noteta-1)*tetap2st
0122           alpha0p2 = -(ptdeltax2+(1.-ptdeltax2)*tplus2)/(1.-2.*tplus2)
0123           alphap2  = (1.+ptdeltax2)*(1.-tplus2)/(1.-2.*tplus2)
0124           betap2   = rkappax2*(1.-tplus2)/(1.-2.*tplus2)**0.5
0125           x(ii)    = rho*(alpha0p2+alphap2*cos(teta))
0126           y(ii)    = -rho*(betap2*sin(teta))
0127           ii       = ii+1
0128         enddo ! noteta
0129 c bottom minus
0130         do noteta=1,nbtetam2
0131           teta     = tetam2in+(noteta-1)*tetam2st
0132           alpha0m2 = -(ptdeltax2-(1.+ptdeltax2)*tmoins2)/(1.-2.*tmoins2)
0133           alpham2  = (1.-ptdeltax2)*(1.-tmoins2)/(1.-2.*tmoins2)
0134           betam2   = rkappax2*(1.-tmoins2)/(1.-2.*tmoins2)**0.5
0135           x(ii)    = rho*(alpha0m2+alpham2*cos(teta))
0136           y(ii)    = -rho*(betam2*sin(teta))
0137 c          if (ii .eq. 150) then
0138            if(iverbose .eq. 1)then
0139               write(msg,*) "y(ii)=",y(ii),betam2,teta,sin(teta),"\n"
0140               call mexPrintf(msg)
0141               endif
0142 
0143 c          endif
0144 
0145           ii       = ii+1
0146         enddo ! noteta
0147 c upper minus
0148         do noteta=1,nbtetam1
0149           teta     = tetam1in+(noteta-1)*tetam1st
0150           alpha0m1 = -(ptdeltax1-(1.+ptdeltax1)*tmoins1)/(1.-2.*tmoins1)
0151           alpham1  = (1.-ptdeltax1)*(1.-tmoins1)/(1.-2.*tmoins1)
0152           betam1   = rkappax1*(1.-tmoins1)/(1.-2.*tmoins1)**0.5
0153           x(ii)    = rho*(alpha0m1+alpham1*cos(teta))
0154           y(ii)    = rho*(betam1*sin(teta))
0155           ii       = ii+1
0156 c      write(7,*)x,y
0157         enddo ! noteta
0158 c upper plus
0159         do noteta=1,nbtetap1
0160           teta     = tetap1in+(noteta-1)*tetap1st
0161           alpha0p1 = -(ptdeltax1+(1.-ptdeltax1)*tplus1)/(1.-2.*tplus1)
0162           alphap1  = (1.+ptdeltax1)*(1.-tplus1)/(1.-2.*tplus1)
0163           betap1   = rkappax1*(1.-tplus1)/(1.-2.*tplus1)**0.5
0164           x(ii)    = rho*(alpha0p1+alphap1*cos(teta))
0165           y(ii)    = -rho*(betap1*sin(teta))
0166           ii       = ii+1
0167 c      write(7,*)x,y
0168         enddo ! noteta
0169       endif ! true/false
0170 
0171 c
0172 c case where the upper part is made of 2 portions of ellipse and the
0173 c bottom part is made of 1 portion of hyperbola (inside) and
0174 c 1 portion of ellipse (outside)
0175       if(mode .eq. 2)then
0176         tetaxp1=0.5*pi-asin(tplus1/(1.-tplus1))
0177         tetap1in=0.
0178         tetap1fi=tetaxp1
0179         nbtetap1=125
0180         tetap1st=(tetap1fi-tetap1in)/(nbtetap1-1)
0181         tetaxm1=0.5*pi+asin(tmoins1/(1.-tmoins1))
0182         tetam1in=tetaxm1
0183         tetam1fi=pi
0184         nbtetam1=125
0185         tetam1st=(tetam1fi-tetam1in)/(nbtetam1-1)
0186         phixm2=dasinh((2.*tmoins2-1.)**0.5/(1.-tmoins2))
0187         if(iverbose .eq. 1)then
0188         write(msg,*) "phixm2=",phixm2,"\n"
0189           call mexPrintf(msg)
0190           endif
0191         phim2in=0.
0192         phim2fi=-phixm2
0193         nbphim2=125
0194         phim2st=(phim2fi-phim2in)/(nbphim2-1)
0195         tetaxp2=0.5*pi-asin(tplus2/(1.-tplus2))
0196         tetap2in=-tetaxp2
0197         tetap2fi=0.
0198         nbtetap2=125
0199         tetap2st=(tetap2fi-tetap2in)/(nbtetap2-1)
0200 c upper plus
0201         ii = 1
0202         do noteta=1,nbtetap1
0203           teta=tetap1in+(noteta-1)*tetap1st
0204           alpha0p1=-(ptdeltax1+(1.-ptdeltax1)*tplus1)/(1.-2.*tplus1)
0205           alphap1=(1.+ptdeltax1)*(1.-tplus1)/(1.-2.*tplus1)
0206           betap1=rkappax1*(1.-tplus1)/(1.-2.*tplus1)**0.5
0207           x(ii) =rho*(alpha0p1+alphap1*cos(teta))
0208           y(ii) =rho*(betap1*sin(teta))
0209           ii = ii+1
0210         enddo ! noteta
0211 c upper minus
0212         do noteta=1,nbtetam1
0213           teta=tetam1in+(noteta-1)*tetam1st
0214           alpha0m1=-(ptdeltax1-(1.+ptdeltax1)*tmoins1)/(1.-2.*tmoins1)
0215           alpham1=(1.-ptdeltax1)*(1.-tmoins1)/(1.-2.*tmoins1)
0216           betam1=rkappax1*(1.-tmoins1)/(1.-2.*tmoins1)**0.5
0217           x(ii)=rho*(alpha0m1+alpham1*cos(teta))
0218           y(ii)=rho*(betam1*sin(teta))
0219           ii = ii+1
0220         enddo ! noteta
0221 c bottom minus
0222         do nophi=1,nbphim2
0223           phi=phim2in+(nophi-1)*phim2st
0224           alpha0m2=-((1.+ptdeltax2)*tmoins2-ptdeltax2)/(2.*tmoins2-1.)
0225           alpham2=(1.-ptdeltax2)*(1.-tmoins2)/(2.*tmoins2-1.)
0226           betam2=rkappax2*(1.-tmoins2)/(2.*tmoins2-1.)**0.5
0227           x(ii)=rho*(alpha0m2+alpham2*cosh(phi))
0228           y(ii)=rho*(betam2*sinh(phi))
0229           ii=ii+1
0230         enddo ! nophi
0231 c bottom plus
0232         do noteta=1,nbtetap2
0233           teta=tetap2in+(noteta-1)*tetap2st
0234           alpha0p2=-(ptdeltax2+(1.-ptdeltax2)*tplus2)/(1.-2.*tplus2)
0235           alphap2=(1.+ptdeltax2)*(1.-tplus2)/(1.-2.*tplus2)
0236           betap2=rkappax2*(1.-tplus2)/(1.-2.*tplus2)**0.5
0237           x(ii)=rho*(alpha0p2+alphap2*cos(teta))
0238           y(ii)=rho*(betap2*sin(teta))
0239           ii = ii+1
0240         enddo ! noteta
0241       endif ! true/false
0242 c
0243 c case where the upper part is made of 1 portion of hyperbola and
0244 c 1 portion of ellipse
0245 c and the bottom part of 2 portions of ellipse
0246 c
0247       if(mode .eq. 3)then
0248          tetaxp2=0.5*pi-asin(tplus2/(1.-tplus2))
0249          tetap2in=0.
0250          tetap2fi=tetaxp2
0251          nbtetap2=125
0252          tetap2st=(tetap2fi-tetap2in)/(nbtetap2-1)
0253          tetaxm2=0.5*pi+asin(tmoins2/(1.-tmoins2))
0254          tetam2in=tetaxm2
0255          tetam2fi=pi
0256          nbtetam2=125
0257          tetam2st=(tetam2fi-tetam2in)/(nbtetam2-1) 
0258          phixm1=dasinh((2.*tmoins1-1.)**0.5/(1.-tmoins1))
0259          if(iverbose .eq. 1)then
0260          write(msg,*) "phixm1=",phixm1,"\n"
0261            call mexPrintf(msg)
0262            endif
0263          phim1in=0.
0264          phim1fi=-phixm1
0265          nbphim1=125
0266          phim1st=(phim1fi-phim1in)/(nbphim1-1)
0267          tetaxp1=0.5*pi-asin(tplus1/(1.-tplus1))
0268          tetap1in=-tetaxp1
0269          tetap1fi=0.
0270          nbtetap1=125
0271          tetap1st=(tetap1fi-tetap1in)/(nbtetap1-1)
0272          ii = 1
0273 c bottom plus
0274          do noteta=1,nbtetap2
0275            teta=tetap2in+(noteta-1)*tetap2st
0276            alpha0p2=-(ptdeltax2+(1.-ptdeltax2)*tplus2)/(1.-2.*tplus2)
0277            alphap2=(1.+ptdeltax2)*(1.-tplus2)/(1.-2.*tplus2)
0278            betap2=rkappax2*(1.-tplus2)/(1.-2.*tplus2)**0.5
0279            x(ii)=rho*(alpha0p2+alphap2*cos(teta))
0280            y(ii)=-rho*(betap2*sin(teta))
0281            ii = ii+1
0282          enddo ! noteta
0283 c bottom minus
0284          do noteta=1,nbtetam2
0285            teta=tetam2in+(noteta-1)*tetam2st
0286            alpha0m2=-(ptdeltax2-(1.+ptdeltax2)*tmoins2)/(1.-2.*tmoins2)
0287            alpham2=(1.-ptdeltax2)*(1.-tmoins2)/(1.-2.*tmoins2)
0288            betam2=rkappax2*(1.-tmoins2)/(1.-2.*tmoins2)**0.5
0289            x(ii)=rho*(alpha0m2+alpham2*cos(teta))
0290            y(ii)=-rho*(betam2*sin(teta))
0291            ii = ii+1
0292          enddo ! noteta
0293 c upper minus
0294          do nophi=1,nbphim1
0295            phi=phim1in+(nophi-1)*phim1st
0296            alpha0m1=-((1.+ptdeltax1)*tmoins1-ptdeltax1)/(2.*tmoins1-1.)
0297            alpham1=(1.-ptdeltax1)*(1.-tmoins1)/(2.*tmoins1-1.)
0298            betam1=rkappax1*(1.-tmoins1)/(2.*tmoins1-1.)**0.5
0299            x(ii)=rho*(alpha0m1+alpham1*cosh(phi))
0300            y(ii)=-rho*(betam1*sinh(phi))
0301            ii = ii+1
0302          enddo ! nophi
0303 c upper plus
0304          do noteta=1,nbtetap1
0305            teta=tetap1in+(noteta-1)*tetap1st
0306            alpha0p1=-(ptdeltax1+(1.-ptdeltax1)*tplus1)/(1.-2.*tplus1)
0307            alphap1=(1.+ptdeltax1)*(1.-tplus1)/(1.-2.*tplus1)
0308            betap1=rkappax1*(1.-tplus1)/(1.-2.*tplus1)**0.5
0309            x(ii)=rho*(alpha0p1+alphap1*cos(teta))
0310            y(ii)=-rho*(betap1*sin(teta))
0311            ii=ii+1
0312          enddo ! noteta
0313       endif ! true/false
0314 c
0315 c case of 2 semi-ellipse with triangularity
0316 c 
0317       if(mode .eq. 4)then
0318         rkappax=0.5*(rkappax1+rkappax2)
0319         ptdeltax=0.5*(ptdeltax1+ptdeltax2)
0320         dzeta0=0.5*(rkappax1-rkappax2)
0321 c external part
0322         tetain=-0.5*pi
0323         tetafi=0.5*pi
0324         nbteta=250
0325         tetast=(tetafi-tetain)/(nbteta-1)
0326         ii = 1
0327         do noteta=1,nbteta
0328           teta=tetain+(noteta-1)*tetast
0329           x(ii)=-ptdeltax+(1.+ptdeltax)*cos(teta)
0330           y(ii)=dzeta0+rkappax*sin(teta)
0331           ii = ii+1
0332         enddo ! noteta
0333 c internal part
0334         tetain=0.5*pi
0335         tetafi=1.5*pi
0336         nbteta=250
0337         tetast=(tetafi-tetain)/(nbteta-1)
0338         do noteta=1,nbteta
0339           teta=tetain+(noteta-1)*tetast
0340           x(ii)=-ptdeltax+(1.-ptdeltax)*cos(teta)
0341           y(ii)=dzeta0+rkappax*sin(teta)
0342           ii = ii +1 
0343         enddo ! noteta
0344       endif ! true/false
0345        sizeOUT = 125*4
0346 C x
0347       if(iverbose .eq. 1)then
0348       write(msg,*) "debut ecriture de x et y\n"
0349       call mexPrintf(msg)
0350       endif
0351       plhs(1) = mxCreateDoubleMatrix(mat500,mat1,0)
0352       x_pr = mxGetPr(plhs(1))
0353       call mxCopyReal8ToPtr(x,x_pr,mat500)
0354 C x
0355       plhs(2) = mxCreateDoubleMatrix(mat500,mat1,0)
0356       y_pr = mxGetPr(plhs(2))
0357       call mxCopyReal8ToPtr(y,y_pr,mat500)
0358 
0359       end
0360 
0361 
0362 
0363 *DECK DASINH
0364       DOUBLE PRECISION FUNCTION DASINH (X)
0365 C***BEGIN PROLOGUE  DASINH
0366 C***PURPOSE  Compute the arc hyperbolic sine.
0367 C***LIBRARY   SLATEC (FNLIB)
0368 C***CATEGORY  C4C
0369 C***TYPE      DOUBLE PRECISION (ASINH-S, DASINH-D, CASINH-C)
0370 C***KEYWORDS  ARC HYPERBOLIC SINE, ASINH, ELEMENTARY FUNCTIONS, FNLIB,
0371 C             INVERSE HYPERBOLIC SINE
0372 C***AUTHOR  Fullerton, W., (LANL)
0373 C***DESCRIPTION
0374 C
0375 C DASINH(X) calculates the double precision arc hyperbolic
0376 C sine for double precision argument X.
0377 C
0378 C***REFERENCES  (NONE)
0379 C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS
0380 C***REVISION HISTORY  (YYMMDD)
0381 C   770601  DATE WRITTEN
0382 C   890531  Changed all specific intrinsics to generic.  (WRB)
0383 C   890531  REVISION DATE from Version 3.2
0384 C   891214  Prologue converted to Version 4.0 format.  (BAB)
0385 C***END PROLOGUE  DASINH
0386       DOUBLE PRECISION X, ASNHCS(39), ALN2, SQEPS, XMAX, Y,
0387      1  DCSEVL, D1MACH
0388       LOGICAL FIRST
0389       SAVE ASNHCS, ALN2, NTERMS, XMAX, SQEPS, FIRST
0390       DATA ASNHCS(  1) / -.1282003991 1738186343 3721273592 68 D+0     /
0391       DATA ASNHCS(  2) / -.5881176118 9951767565 2117571383 62 D-1     /
0392       DATA ASNHCS(  3) / +.4727465432 2124815640 7252497560 29 D-2     /
0393       DATA ASNHCS(  4) / -.4938363162 6536172101 3601747902 73 D-3     /
0394       DATA ASNHCS(  5) / +.5850620705 8557412287 4948352593 21 D-4     /
0395       DATA ASNHCS(  6) / -.7466998328 9313681354 7550692171 88 D-5     /
0396       DATA ASNHCS(  7) / +.1001169358 3558199265 9661920158 12 D-5     /
0397       DATA ASNHCS(  8) / -.1390354385 8708333608 6164722588 86 D-6     /
0398       DATA ASNHCS(  9) / +.1982316948 3172793547 3173602371 48 D-7     /
0399       DATA ASNHCS( 10) / -.2884746841 7848843612 7472728003 17 D-8     /
0400       DATA ASNHCS( 11) / +.4267296546 7159937953 4575149959 07 D-9     /
0401       DATA ASNHCS( 12) / -.6397608465 4366357868 7526323096 81 D-10    /
0402       DATA ASNHCS( 13) / +.9699168608 9064704147 8782931311 79 D-11    /
0403       DATA ASNHCS( 14) / -.1484427697 2043770830 2466583656 96 D-11    /
0404       DATA ASNHCS( 15) / +.2290373793 9027447988 0401843789 83 D-12    /
0405       DATA ASNHCS( 16) / -.3558839513 2732645159 9789426513 10 D-13    /
0406       DATA ASNHCS( 17) / +.5563969408 0056789953 3745390885 54 D-14    /
0407       DATA ASNHCS( 18) / -.8746250959 9624678045 6665935201 62 D-15    /
0408       DATA ASNHCS( 19) / +.1381524884 4526692155 8688022981 29 D-15    /
0409       DATA ASNHCS( 20) / -.2191668828 2900363984 9551422641 49 D-16    /
0410       DATA ASNHCS( 21) / +.3490465852 4827565638 3139237068 80 D-17    /
0411       DATA ASNHCS( 22) / -.5578578840 0895742439 6301570321 06 D-18    /
0412       DATA ASNHCS( 23) / +.8944514661 7134012551 0508827989 33 D-19    /
0413       DATA ASNHCS( 24) / -.1438342634 6571317305 5518452394 66 D-19    /
0414       DATA ASNHCS( 25) / +.2319181187 2169963036 3261446826 66 D-20    /
0415       DATA ASNHCS( 26) / -.3748700795 3314343674 5706045439 99 D-21    /
0416       DATA ASNHCS( 27) / +.6073210982 2064279404 5492428800 00 D-22    /
0417       DATA ASNHCS( 28) / -.9859940276 4633583177 3701734400 00 D-23    /
0418       DATA ASNHCS( 29) / +.1603921745 2788496315 2326382933 33 D-23    /
0419       DATA ASNHCS( 30) / -.2613884735 0287686596 7161343999 99 D-24    /
0420       DATA ASNHCS( 31) / +.4267084960 6857390833 3581653333 33 D-25    /
0421       DATA ASNHCS( 32) / -.6977021703 9185243299 7307733333 33 D-26    /
0422       DATA ASNHCS( 33) / +.1142508833 6806858659 8126933333 33 D-26    /
0423       DATA ASNHCS( 34) / -.1873529207 8860968933 0210133333 33 D-27    /
0424       DATA ASNHCS( 35) / +.3076358441 4464922794 0659200000 00 D-28    /
0425       DATA ASNHCS( 36) / -.5057736403 1639824787 0463999999 99 D-29    /
0426       DATA ASNHCS( 37) / +.8325075471 2689142224 2133333333 33 D-30    /
0427       DATA ASNHCS( 38) / -.1371845728 2501044163 9253333333 33 D-30    /
0428       DATA ASNHCS( 39) / +.2262986842 6552784104 1066666666 66 D-31    /
0429       DATA ALN2 / 0.6931471805 5994530941 7232121458 18D0 /
0430       DATA FIRST /.TRUE./
0431 C***FIRST EXECUTABLE STATEMENT  DASINH
0432       IF (FIRST) THEN
0433          NTERMS = INITDS (ASNHCS, 39, 0.1*REAL(D1MACH(3)) )
0434          SQEPS = SQRT(D1MACH(3))
0435          XMAX = 1.0D0/SQEPS
0436       ENDIF
0437       FIRST = .FALSE.
0438 C
0439       Y = ABS(X)
0440       IF (Y.GT.1.0D0) GO TO 20
0441 C
0442       DASINH = X
0443       IF (Y.GT.SQEPS) DASINH = X*(1.0D0 + DCSEVL (2.D0*X*X-1.D0,
0444      1  ASNHCS, NTERMS) )
0445       RETURN
0446  20   IF (Y.LT.XMAX) DASINH = LOG (Y+SQRT(Y*Y+1.D0))
0447       IF (Y.GE.XMAX) DASINH = ALN2 + LOG(Y)
0448       DASINH = SIGN (DASINH, X)
0449       RETURN
0450 C
0451       END
0452 
0453 
0454             DOUBLE PRECISION FUNCTION D1MACH(I)
0455       integer I
0456 C
0457 C  DOUBLE-PRECISION MACHINE CONSTANTS
0458 C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
0459 C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
0460 C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
0461 C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
0462 C  D1MACH( 5) = LOG10(B)
0463 C
0464       integer SMALL(2)
0465       integer LARGE(2)
0466       integer RIGHT(2)
0467       integer DIVER(2)
0468       integer LOG10(2)
0469       integer SC, CRAY1(38), J
0470       COMMON /D9MACH/ CRAY1
0471       SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC
0472       DOUBLE PRECISION DMACH(5)
0473       EQUIVALENCE (DMACH(1),SMALL(1))
0474       EQUIVALENCE (DMACH(2),LARGE(1))
0475       EQUIVALENCE (DMACH(3),RIGHT(1))
0476       EQUIVALENCE (DMACH(4),DIVER(1))
0477       EQUIVALENCE (DMACH(5),LOG10(1))
0478 C  THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES.
0479 C  R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF
0480 C  D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR
0481 C  MANY MACHINES YET.
0482 C  TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1
0483 C  ON THE NEXT LINE
0484       DATA SC/0/
0485 C  AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW.
0486 C  CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY
0487 C          mail netlib@research.bell-labs.com
0488 C          send old1mach from blas
0489 C  PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com.
0490 C
0491 C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
0492 C      DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
0493 C      DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
0494 C      DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
0495 C      DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
0496 C      DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/
0497 C
0498 C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
0499 C     32-BIT INTEGERS.
0500 C      DATA SMALL(1),SMALL(2) /    8388608,           0 /
0501 C      DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
0502 C      DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
0503 C      DATA DIVER(1),DIVER(2) /  620756992,           0 /
0504 C      DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/
0505 C
0506 C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
0507 C      DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
0508 C      DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
0509 C      DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
0510 C      DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
0511 C      DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/
0512 C
0513 C     ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES.
0514       IF (SC .NE. 987) THEN
0515          DMACH(1) = 1.D13
0516          IF (      SMALL(1) .EQ. 1117925532
0517      *       .AND. SMALL(2) .EQ. -448790528) THEN
0518 *           *** IEEE BIG ENDIAN ***
0519             SMALL(1) = 1048576
0520             SMALL(2) = 0
0521             LARGE(1) = 2146435071
0522             LARGE(2) = -1
0523             RIGHT(1) = 1017118720
0524             RIGHT(2) = 0
0525             DIVER(1) = 1018167296
0526             DIVER(2) = 0
0527             LOG10(1) = 1070810131
0528             LOG10(2) = 1352628735
0529          ELSE IF ( SMALL(2) .EQ. 1117925532
0530      *       .AND. SMALL(1) .EQ. -448790528) THEN
0531 *           *** IEEE LITTLE ENDIAN ***
0532             SMALL(2) = 1048576
0533             SMALL(1) = 0
0534             LARGE(2) = 2146435071
0535             LARGE(1) = -1
0536             RIGHT(2) = 1017118720
0537             RIGHT(1) = 0
0538             DIVER(2) = 1018167296
0539             DIVER(1) = 0
0540             LOG10(2) = 1070810131
0541             LOG10(1) = 1352628735
0542          ELSE IF ( SMALL(1) .EQ. -2065213935
0543      *       .AND. SMALL(2) .EQ. 10752) THEN
0544 *               *** VAX WITH D_FLOATING ***
0545             SMALL(1) = 128
0546             SMALL(2) = 0
0547             LARGE(1) = -32769
0548             LARGE(2) = -1
0549             RIGHT(1) = 9344
0550             RIGHT(2) = 0
0551             DIVER(1) = 9472
0552             DIVER(2) = 0
0553             LOG10(1) = 546979738
0554             LOG10(2) = -805796613
0555          ELSE IF ( SMALL(1) .EQ. 1267827943
0556      *       .AND. SMALL(2) .EQ. 704643072) THEN
0557 *               *** IBM MAINFRAME ***
0558             SMALL(1) = 1048576
0559             SMALL(2) = 0
0560             LARGE(1) = 2147483647
0561             LARGE(2) = -1
0562             RIGHT(1) = 856686592
0563             RIGHT(2) = 0
0564             DIVER(1) = 873463808
0565             DIVER(2) = 0
0566             LOG10(1) = 1091781651
0567             LOG10(2) = 1352628735
0568          ELSE IF ( SMALL(1) .EQ. 1120022684
0569      *       .AND. SMALL(2) .EQ. -448790528) THEN
0570 *           *** CONVEX C-1 ***
0571             SMALL(1) = 1048576
0572             SMALL(2) = 0
0573             LARGE(1) = 2147483647
0574             LARGE(2) = -1
0575             RIGHT(1) = 1019215872
0576             RIGHT(2) = 0
0577             DIVER(1) = 1020264448
0578             DIVER(2) = 0
0579             LOG10(1) = 1072907283
0580             LOG10(2) = 1352628735
0581          ELSE IF ( SMALL(1) .EQ. 815547074
0582      *       .AND. SMALL(2) .EQ. 58688) THEN
0583 *           *** VAX G-FLOATING ***
0584             SMALL(1) = 16
0585             SMALL(2) = 0
0586             LARGE(1) = -32769
0587             LARGE(2) = -1
0588             RIGHT(1) = 15552
0589             RIGHT(2) = 0
0590             DIVER(1) = 15568
0591             DIVER(2) = 0
0592             LOG10(1) = 1142112243
0593             LOG10(2) = 2046775455
0594          ELSE
0595             DMACH(2) = 1.D27 + 1
0596             DMACH(3) = 1.D27
0597             LARGE(2) = LARGE(2) - RIGHT(2)
0598             IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN
0599                CRAY1(1) = 67291416
0600                DO 10 J = 1, 20
0601                   CRAY1(J+1) = CRAY1(J) + CRAY1(J)
0602  10               CONTINUE
0603                CRAY1(22) = CRAY1(21) + 321322
0604                DO 20 J = 22, 37
0605                   CRAY1(J+1) = CRAY1(J) + CRAY1(J)
0606  20               CONTINUE
0607                IF (CRAY1(38) .EQ. SMALL(1)) THEN
0608 *                  *** CRAY ***
0609                   CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0)
0610                   SMALL(2) = 0
0611                   CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215)
0612                   CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214)
0613                   CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0)
0614                   RIGHT(2) = 0
0615                   CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0)
0616                   DIVER(2) = 0
0617                   CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215)
0618                   CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388)
0619                ELSE
0620                   WRITE(*,9000)
0621                   STOP 779
0622                   END IF
0623             ELSE
0624                WRITE(*,9000)
0625                STOP 779
0626                END IF
0627             END IF
0628          SC = 987
0629          END IF
0630 *    SANITY CHECK
0631       IF (DMACH(4) .GE. 1.0D0) STOP 778
0632       IF (I .LT. 1 .OR. I .GT. 5) THEN
0633       if(iverbose .eq. 1)then
0634       write(msg,*) "D1MACH(I): I =",I," is out of bounds.\n"
0635       call mexPrintf(msg)
0636       endif
0637          STOP
0638          END IF
0639       D1MACH = DMACH(I)
0640       RETURN
0641  9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/
0642      *' appropriate for your machine.')
0643 * /* Standard C source for D1MACH -- remove the * in column 1 */
0644 *#include 
0645 *#include 
0646 *#include 
0647 *double d1mach_(long *i)
0648 *{
0649 *  switch(*i){
0650 *    case 1: return DBL_MIN;
0651 *    case 2: return DBL_MAX;
0652 *    case 3: return DBL_EPSILON/FLT_RADIX;
0653 *    case 4: return DBL_EPSILON;
0654 *    case 5: return log10((double)FLT_RADIX);
0655 *    }
0656 *  fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
0657 *  exit(1); return 0; /* some compilers demand return values */
0658 *}
0659       END
0660       SUBROUTINE I1MCRY(A, A1, B, C, D)
0661 **** SPECIAL COMPUTATION FOR OLD CRAY MACHINES ****
0662       integer A, A1, B, C, D
0663       A1 = 16777216*B + C
0664       A = 16777216*A1 + D
0665       END
0666 
0667 
0668 
0669 *DECK DCSEVL
0670       DOUBLE PRECISION FUNCTION DCSEVL (X, CS, N)
0671 C***BEGIN PROLOGUE  DCSEVL
0672 C***PURPOSE  Evaluate a Chebyshev series.
0673 C***LIBRARY   SLATEC (FNLIB)
0674 C***CATEGORY  C3A2
0675 C***TYPE      DOUBLE PRECISION (CSEVL-S, DCSEVL-D)
0676 C***KEYWORDS  CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS
0677 C***AUTHOR  Fullerton, W., (LANL)
0678 C***DESCRIPTION
0679 C
0680 C  Evaluate the N-term Chebyshev series CS at X.  Adapted from
0681 C  a method presented in the paper by Broucke referenced below.
0682 C
0683 C       Input Arguments --
0684 C  X    value at which the series is to be evaluated.
0685 C  CS   array of N terms of a Chebyshev series.  In evaluating
0686 C       CS, only half the first coefficient is summed.
0687 C  N    number of terms in array CS.
0688 C
0689 C***REFERENCES  R. Broucke, Ten subroutines for the manipulation of
0690 C                 Chebyshev series, Algorithm 446, Communications of
0691 C                 the A.C.M. 16, (1973) pp. 254-256.
0692 C               L. Fox and I. B. Parker, Chebyshev Polynomials in
0693 C                 Numerical Analysis, Oxford University Press, 1968,
0694 C                 page 56.
0695 C***ROUTINES CALLED  D1MACH, XERMSG
0696 C***REVISION HISTORY  (YYMMDD)
0697 C   770401  DATE WRITTEN
0698 C   890831  Modified array declarations.  (WRB)
0699 C   890831  REVISION DATE from Version 3.2
0700 C   891214  Prologue converted to Version 4.0 format.  (BAB)
0701 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
0702 C   900329  Prologued revised extensively and code rewritten to allow
0703 C           X to be slightly outside interval (-1,+1).  (WRB)
0704 C   920501  Reformatted the REFERENCES section.  (WRB)
0705 C***END PROLOGUE  DCSEVL
0706       DOUBLE PRECISION B0, B1, B2, CS(*), ONEPL, TWOX, X, D1MACH
0707       LOGICAL FIRST
0708       SAVE FIRST, ONEPL
0709       DATA FIRST /.TRUE./
0710 C***FIRST EXECUTABLE STATEMENT  DCSEVL
0711       IF (FIRST) ONEPL = 1.0D0 + D1MACH(4)
0712       FIRST = .FALSE.
0713 c      IF (N .LT. 1) PRINT *,'SLATEC', 'DCSEVL','SLATEC', 'DCSEVL',
0714 c     +   'NUMBER OF TERMS .LE. 0', 2, 2
0715 c      IF (N .GT. 1000) PRINT *,'SLATEC', 'DCSEVL','SLATEC', 'DCSEVL',
0716 c     +   'NUMBER OF TERMS .GT. 1000', 3, 2
0717 c      IF (ABS(X) .GT. ONEPL) PRINT *,'SLATEC', 'DCSEVL','SLATEC', 'DCSEVL',
0718 c     +   'X OUTSIDE THE INTERVAL (-1,+1)', 1, 1
0719 C
0720       B1 = 0.0D0
0721       B0 = 0.0D0
0722       TWOX = 2.0D0*X
0723       DO 10 I = 1,N
0724          B2 = B1
0725          B1 = B0
0726          NI = N + 1 - I
0727          B0 = TWOX*B1 - B2 + CS(NI)
0728    10 CONTINUE
0729 C
0730       DCSEVL = 0.5D0*(B0-B2)
0731 C
0732       RETURN
0733       END
0734 
0735 
0736 *DECK INITDS
0737       FUNCTION INITDS (OS, NOS, ETA)
0738 C***BEGIN PROLOGUE  INITDS
0739 C***PURPOSE  Determine the number of terms needed in an orthogonal
0740 C            polynomial series so that it meets a specified accuracy.
0741 C***LIBRARY   SLATEC (FNLIB)
0742 C***CATEGORY  C3A2
0743 C***TYPE      DOUBLE PRECISION (INITS-S, INITDS-D)
0744 C***KEYWORDS  CHEBYSHEV, FNLIB, INITIALIZE, ORTHOGONAL POLYNOMIAL,
0745 C             ORTHOGONAL SERIES, SPECIAL FUNCTIONS
0746 C***AUTHOR  Fullerton, W., (LANL)
0747 C***DESCRIPTION
0748 C
0749 C  Initialize the orthogonal series, represented by the array OS, so
0750 C  that INITDS is the number of terms needed to insure the error is no
0751 C  larger than ETA.  Ordinarily, ETA will be chosen to be one-tenth
0752 C  machine precision.
0753 C
0754 C             Input Arguments --
0755 C   OS     double precision array of NOS coefficients in an orthogonal
0756 C          series.
0757 C   NOS    number of coefficients in OS.
0758 C   ETA    single precision scalar containing requested accuracy of
0759 C          series.
0760 C
0761 C***REFERENCES  (NONE)
0762 C***ROUTINES CALLED  XERMSG
0763 C***REVISION HISTORY  (YYMMDD)
0764 C   770601  DATE WRITTEN
0765 C   890531  Changed all specific intrinsics to generic.  (WRB)
0766 C   890831  Modified array declarations.  (WRB)
0767 C   891115  Modified error message.  (WRB)
0768 C   891115  REVISION DATE from Version 3.2
0769 C   891214  Prologue converted to Version 4.0 format.  (BAB)
0770 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
0771 C***END PROLOGUE  INITDS
0772       DOUBLE PRECISION OS(*)
0773 C***FIRST EXECUTABLE STATEMENT  INITDS
0774 c      IF (NOS .LT. 1) PRINT *,'SLATEC', 'INITDS','SLATEC', 'INITDS',
0775 c     +   'Number of coefficients is less than 1', 2, 1
0776 C
0777       ERR = 0.
0778       DO 10 II = 1,NOS
0779         I = NOS + 1 - II
0780         ERR = ERR + ABS(REAL(OS(I)))
0781         IF (ERR.GT.ETA) GO TO 20
0782    10 CONTINUE
0783 C
0784    20  continue
0785 c      IF (I .EQ. NOS) PRINT *,'SLATEC', 'INITDS','SLATEC', 'INITDS',
0786 c     +   'Chebyshev series too short for specified accuracy', 1, 1
0787       INITDS = I
0788 C
0789       RETURN
0790       END
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.