c ************************************************************ c c chapman function for grazing incidence c c ************************************************************ function ch(x,z) integer ierr real ch,x double precision ch1,sublim,z,error,aerr,rerr double precision dcadre,fch external fch,dcadre common/integ/ax,az double precision ax,az c ax=x az=z aerr= 0.0d0 rerr= 0.0001d0 sublim = 1.0d-4 ch1 = dcadre( fch, sublim, az, aerr, rerr, error, ierr ) ch = ax * sin(az) * ch1 return end c ************************************************************ c c integrand for the chapman function c c ************************************************************ double precision function fch(gi) double precision gi common/integ/x,z double precision x,z ! real x,z,a fch = (1./sin(gi))**2. * exp( x *(1.- sin(z)/sin(gi)) ) return end ****************************************************************************** double precision function dcadre (f,a,b,aerr,rerr,error,ier) c specifications for arguments integer ier double precision f,a,b,aerr,rerr,error c specifications for local variables integer ibegs(30),maxts,maxtbl,mxstge,ibeg,ii,nnleft integer i,n2,iii,istep2,iend,istep,l,lm1,it,istage,n double precision t(10,10),r(10),ait(10),dif(10),rn(4),ts(2049) double precision begin(30),finis(30),est(30) double precision h2tol,aittol,length,jumptl,zero,p1,half,one double precision two,four,fourp5,ten,hun,cadre,aitlow double precision stepmn,stepnm,stage,curest,fnsize,hrerr double precision prever,beg,fbeg,edn,fend,step,astep,tabs,hovn double precision fn,sum,sumabs,vint,tabtlm,ergl,ergoal double precision erra,errr,fextrp,errer,diff,sing,fextm1 double precision h2next,singnx,slope,fbeg2,erret,h2tfex,fi logical h2conv,aitken,right,reglar,reglsv(30) data aitlow,h2tol,aittol,jumptl,maxts,maxtbl,mxstge 1 /1.1d0,.15d0,.1d0,.01d0,2049,10,30/ data rn(1),rn(2),rn(3),rn(4)/ 1 .7142005d0,.3466282d0,.843751d0,.1263305d0/ data zero,p1,half,one,two,four,fourp5,ten,hun 1 /0.0d0,0.1d0,0.5d0,1.0d0,2.0d0,4.0d0, 2 4.5d0,10.0d0,100.0d0/ c first executable statement ier = 0 cadre = zero error = zero curest = zero vint = zero length = dabs(b-a) if (length .eq. zero) go to 215 if (rerr .gt. p1 .or. rerr .lt. zero) go to 210 hrerr = rerr+hun if (aerr .eq. zero .and. hrerr .le. hun) go to 210 errr = rerr erra = dabs(aerr) stepmn = length/(two**mxstge) stepnm = dmax1(length,dabs(a),dabs(b))*ten stage = half istage = 1 fnsize = zero prever = zero reglar = .false. c the given interval of integration c is the first interval considered. beg = a fbeg = f(beg)*half ts(1) = fbeg ibeg = 1 edn = b fend = f(edn)*half ts(2) = fend iend = 2 5 right = .false. c investigation of a particular c subinterval begins at this point. 10 step = edn - beg astep = dabs(step) if (astep .lt. stepmn) go to 205 hrerr = stepnm+astep if (hrerr .eq. stepnm) go to 205 t(1,1) = fbeg + fend tabs = dabs(fbeg) + dabs(fend) l = 1 n = 1 h2conv = .false. aitken = .false. 15 lm1 = l l = l + 1 c calculate the next trapezoid sum, c t(l,1), which is based on *n2* + 1 c equispaced points. here, c n2 = n*2 = 2**(l-1). n2 = n+n fn = n2 istep = (iend - ibeg)/n if (istep .gt. 1) go to 25 ii = iend iend = iend + n if (iend .gt. maxts) go to 200 hovn = step/fn iii = iend fi = one do 20 i=1,n2,2 ts(iii) = ts(ii) ts(iii-1) = f(edn - fi * hovn) fi = fi+two iii = iii-2 ii = ii-1 20 continue istep = 2 25 istep2 = ibeg + istep/2 sum = zero sumabs = zero do 30 i=istep2,iend,istep sum = sum + ts(i) sumabs = sumabs + dabs(ts(i)) 30 continue t(l,1) = t(l-1,1)*half+sum/fn tabs = tabs*half+sumabs/fn n = n2 c get preliminary value for *vint* c from last trapezoid sum and update c the error requirement *ergoal* c for this subinterval. it = 1 vint = step*t(l,1) tabtlm = tabs*ten fnsize = dmax1(fnsize,dabs(t(l,1))) ergl = astep*fnsize*ten ergoal = stage*dmax1(erra,errr*dabs(curest+vint)) c complete row l and column l of *t* c array. fextrp = one do 35 i=1,lm1 fextrp = fextrp*four t(i,l) = t(l,i) - t(l-1,i) t(l,i+1) = t(l,i) + t(i,l)/(fextrp-one) 35 continue errer = astep*dabs(t(1,l)) c preliminary decision procedure c if l = 2 and t(2,1) = t(1,1), c go to 135 to follow up the c impression that intergrand is c straight line. if (l .gt. 2) go to 40 hrerr = tabs+p1*dabs(t(1,2)) if (hrerr .eq. tabs) go to 135 go to 15 c caculate next ratios for c columns 1,...,l-2 of t-table c ratio is set to zero if difference c in last two entries of column is c about zero 40 do 45 i=2,lm1 diff = zero hrerr = tabtlm+dabs(t(i-1,l)) if (hrerr .ne. tabtlm) diff = t(i-1,lm1)/t(i-1,l) t(i-1,lm1) = diff 45 continue if (dabs(four-t(1,lm1)) .le. h2tol) go to 60 if (t(1,lm1) .eq. zero) go to 55 if (dabs(two-dabs(t(1,lm1))) .lt. jumptl) go to 130 if (l .eq. 3) go to 15 h2conv = .false. if (dabs((t(1,lm1)-t(1,l-2))/t(1,lm1)) .le. aittol) go to 75 50 if (reglar) go to 55 if (l .eq. 4) go to 15 hrerr = ergl+errer 55 if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 175 go to 145 c cautious romberg extrapolation 60 if (h2conv) go to 65 aitken = .false. h2conv = .true. 65 fextrp = four 70 it = it + 1 vint = step*t(l,it) errer = dabs(step/(fextrp-one)*t(it-1,l)) if (errer .le. ergoal) go to 160 hrerr = ergl+errer if (hrerr .eq. ergl) go to 160 if (it .eq. lm1) go to 125 if (t(it,lm1) .eq. zero) go to 70 if (t(it,lm1) .le. fextrp) go to 125 if (dabs(t(it,lm1)/four-fextrp)/fextrp .lt. aittol) 1 fextrp = fextrp*four go to 70 c integrand may have x**alpha type c singularity c resulting in a ratio of *sing* = c 2**(alpha + 1) 75 if (t(1,lm1) .lt. aitlow) go to 175 if (aitken) go to 80 h2conv = .false. aitken = .true. 80 fextrp = t(l-2,lm1) if (fextrp .gt. fourp5) go to 65 if (fextrp .lt. aitlow) go to 175 if (dabs(fextrp-t(l-3,lm1))/t(1,lm1) .gt. h2tol) go to 175 sing = fextrp fextm1 = one/(fextrp - one) ait(1) = zero do 85 i=2,l ait(i) = t(i,1) + (t(i,1)-t(i-1,1))*fextm1 r(i) = t(1,i-1) dif(i) = ait(i) - ait(i-1) 85 continue it = 2 90 vint = step*ait(l) errer = errer*fextm1 hrerr = ergl+errer if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 95 ier = max0(ier,65) go to 160 95 it = it + 1 if (it .eq. lm1) go to 125 if (it .gt. 3) go to 100 h2next = four singnx = sing+sing 100 if (h2next .lt. singnx) go to 105 fextrp = singnx singnx = singnx+singnx go to 110 105 fextrp = h2next h2next = four*h2next 110 do 115 i=it,lm1 r(i+1) = zero hrerr = tabtlm+dabs(dif(i+1)) if (hrerr .ne. tabtlm) r(i+1) = dif(i)/dif(i+1) 115 continue h2tfex = -h2tol*fextrp if (r(l) - fextrp .lt. h2tfex) go to 125 if (r(l-1)-fextrp .lt. h2tfex) go to 125 errer = astep*dabs(dif(l)) fextm1 = one/(fextrp - one) do 120 i=it,l ait(i) = ait(i) + dif(i)*fextm1 dif(i) = ait(i) - ait(i-1) 120 continue go to 90 c current trapezoid sum and resulting c extrapolated values did not give c a small enough *errer*. c note -- having prever .lt. errer c is an almost certain sign of c beginning trouble with in the func- c tion values. hence, a watch for, c and control of, noise should c begin here. 125 fextrp = dmax1(prever/errer,aitlow) prever = errer if (l .lt. 5) go to 15 if (l-it .gt. 2 .and. istage .lt. mxstge) go to 170 erret = errer/(fextrp**(maxtbl-l)) hrerr = ergl+erret if (erret .gt. ergoal .and. hrerr .ne. ergl) go to 170 go to 15 c integrand has jump (see notes) 130 hrerr = ergl+errer if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 170 c note that 2*fn = 2**l diff = dabs(t(1,l))*(fn+fn) go to 160 c integrand is straight line c test this assumption by comparing c the value of the integrand at c four *randomly chosen* points with c the value of the straight line c interpolating the integrand at the c two end points of the sub-interval. c if test is passed, accept *vint* 135 slope = (fend-fbeg)*two fbeg2 = fbeg+fbeg do 140 i=1,4 diff = dabs(f(beg+rn(i)*step) - fbeg2-rn(i)*slope) hrerr = tabtlm+diff if(hrerr .ne. tabtlm) go to 155 140 continue go to 160 c noise may be dominant feature c estimate noise level by comparing c the value of the integrand at c four *randomly chosen* points with c the value of the straight line c interpolating the integrand at the c two endpoints. if small enough, c accept *vint* 145 slope = (fend-fbeg)*two fbeg2 = fbeg+fbeg i = 1 150 diff = dabs(f(beg+rn(i)*step) - fbeg2-rn(i)*slope) 155 errer = dmax1(errer,astep*diff) hrerr = ergl+errer if (errer .gt. ergoal .and. hrerr .ne. ergl) go to 175 i = i+1 if (i .le. 4) go to 150 ier = 66 c intergration over current sub- c interval successful c add *vint* to *dcadre* and *errer* c to *error*, then set up next sub- c interval, if any. 160 cadre = cadre + vint error = error + errer if (right) go to 165 istage = istage - 1 if (istage .eq. 0) go to 220 reglar = reglsv(istage) beg = begin(istage) edn = finis(istage) curest = curest - est(istage+1) + vint iend = ibeg - 1 fend = ts(iend) ibeg = ibegs(istage) go to 180 165 curest = curest + vint stage = stage+stage iend = ibeg ibeg = ibegs(istage) edn = beg beg = begin(istage) fend = fbeg fbeg = ts(ibeg) go to 5 c integration over current subinterval c is unsuccessful. mark subinterval c for further subdivision. set up c next subinterval. 170 reglar = .true. 175 if (istage .eq. mxstge) go to 205 if (right) go to 185 reglsv(istage+1) = reglar begin(istage) = beg ibegs(istage) = ibeg stage = stage*half 180 right = .true. beg = (beg+edn)*half ibeg = (ibeg+iend)/2 ts(ibeg) = ts(ibeg)*half fbeg = ts(ibeg) go to 10 185 nnleft = ibeg - ibegs(istage) if (iend+nnleft .ge. maxts) go to 200 iii = ibegs(istage) ii = iend do 190 i=iii,ibeg ii = ii + 1 ts(ii) = ts(i) 190 continue do 195 i=ibeg,ii ts(iii) = ts(i) iii = iii + 1 195 continue iend = iend + 1 ibeg = iend - nnleft fend = fbeg fbeg = ts(ibeg) finis(istage) = edn edn = beg beg = begin(istage) begin(istage) = edn reglsv(istage) = reglar istage = istage + 1 reglar = reglsv(istage) est(istage) = vint curest = curest + est(istage) go to 5 c failure to handle given integra- c tion problem 200 ier = 131 go to 215 205 ier = 132 go to 215 210 ier = 133 215 cadre = curest + vint 220 dcadre = cadre 9000 continue if (ier .ne. 0) call uertst (ier,6hdcadre) 9005 return end ****************************************************************************** subroutine uertst (ier,name) c specifications for arguments integer ier integer name(2) c specifications for local variables integer i,ieq,ieqdf,iounit,level,levold,nameq(6), * namset(6),namupk(6),nin,nmtb data namset/1hu,1he,1hr,1hs,1he,1ht/ data nameq/6*1h / data level/4/,ieqdf/0/,ieq/1h=/ c unpack name into namupk c first executable statement call uspkd (name,6,namupk,nmtb) c get output unit number call ugetio(1,nin,iounit) c check ier if (ier.gt.999) go to 25 if (ier.lt.-32) go to 55 if (ier.le.128) go to 5 if (level.lt.1) go to 30 c print terminal message if (ieqdf.eq.1) write(iounit,35) ier,nameq,ieq,namupk if (ieqdf.eq.0) write(iounit,35) ier,namupk go to 30 5 if (ier.le.64) go to 10 if (level.lt.2) go to 30 c print warning with fix message c if (ieqdf.eq.1) write(iounit,40) ier,nameq,ieq,namupk c if (ieqdf.eq.0) write(iounit,40) ier,namupk if (ieqdf.eq.1) continue if (ieqdf.eq.0) continue go to 30 10 if (ier.le.32) go to 15 c print warning message if (level.lt.3) go to 30 if (ieqdf.eq.1) write(iounit,45) ier,nameq,ieq,namupk if (ieqdf.eq.0) write(iounit,45) ier,namupk go to 30 15 continue c check for uerset call do 20 i=1,6 if (namupk(i).ne.namset(i)) go to 25 20 continue levold = level level = ier ier = levold if (level.lt.0) level = 4 if (level.gt.4) level = 4 go to 30 25 continue if (level.lt.4) go to 30 c print non-defined message if (ieqdf.eq.1) write(iounit,50) ier,nameq,ieq,namupk if (ieqdf.eq.0) write(iounit,50) ier,namupk 30 ieqdf = 0 return 35 format(19h *** terminal error,10x,7h(ier = ,i3, 1 20h) from imsl routine ,6a1,a1,6a1) 40 format(27h *** warning with fix error,2x,7h(ier = ,i3, 1 20h) from imsl routine ,6a1,a1,6a1) 45 format(18h *** warning error,11x,7h(ier = ,i3, 1 20h) from imsl routine ,6a1,a1,6a1) 50 format(20h *** undefined error,9x,7h(ier = ,i5, 1 20h) from imsl routine ,6a1,a1,6a1) c c save p for p = r case c p is the page namupk c r is the routine namupk 55 ieqdf = 1 do 60 i=1,6 60 nameq(i) = namupk(i) 65 return end ************************************************************************ subroutine uspkd (packed,nchars,unpakd,nchmtb) c specifications for arguments integer nc,nchars,nchmtb c c integer unpakd(1),iblank integer unpakd(nchars),iblank character packed(1) ! Modificado el 29-Ago-2001 porque ! integer packed(1) esto deberia ser CHARACTER, no?? data iblank /1h / c initialize nchmtb nchmtb = 0 c return if nchars is le zero if(nchars.le.0) return c set nc=number of chars to be decoded nc = min0 (129,nchars) ! decode (nc,150,packed) (unpakd(i),i=1,nc) read (packed,150) (unpakd(i),i=1,nc) 150 format (129a1) c check unpakd array and set nchmtb c based on trailing blanks found do 200 n = 1,nc nn = nc - n + 1 if(unpakd(nn) .ne. 536870912) go to 210 200 continue 210 nchmtb = nn return end **************************************************************************** subroutine ugetio(iopt,nin,nout) c specifications for arguments integer iopt,nin,nout c specifications for local variables integer nind,noutd data nind/5/,noutd/6/ c first executable statement if (iopt.eq.3) go to 10 if (iopt.eq.2) go to 5 if (iopt.ne.1) go to 9005 nin = nind nout = noutd go to 9005 5 nind = nin go to 9005 10 noutd = nout 9005 return end