c******************************************************************* c Program solves Possio's integral equation for the pressure c distribution and the aerodynamic coefficients of the pitching c and plunging motion in subsonic flow c c Programming and Copyright: V. Carstens (DLR-AE) c c V. Carstens, Berechnung der instationaeren Druck- c verteilung an harmonisch schwingenden Gittern in ebener c Unterschallstroemung, Teil II, DFVLR IB 253-75J02 c DFVLR-AVA, Goettingen 1975 c c Last revision: 15.09.2001 c c Modified for use in FLT2DQ with kind permission of c the author by W. Send (DLR-AE) - 17.01.2005 c Some parts of the program are not in use in this application c******************************************************************* c Subroutine POSSIOF implicit real*8(a-b,d-h,o-z),complex*16(c) real*8 mach,machq,mue character*20 outdat Common /POSSIO/cqval,p_mach,p_redf,p_xr,p_xrud,ifvek,idr,outdat c parameter(nx=51,nmatx=10*nx+1) c dimension cmat(nx,nx),ckof(nx),cw(nx) dimension cvek(nmatx),fmjac(nx,nmatx),gew(nmatx) dimension ifvek(3),CQVAL(3,3) c common /piwert/ pih,pi common /imag/ ci common /ma/ mach,machq,beta,beta2 common /rf/ redf common /gam/ mue,gamma,alpha common /euz/ euzahl common /uint/ cu0 c save c c************************************************************ c Constants c************************************************************ c pih = asin(1.d0) pi = 2.d0*pih euzahl = 0.5772156649015328d0 ci = (0.d0,1.d0) c c************************************************************ c Reading the input data replace by common /POSSIO/ c************************************************************ c mach = p_mach redf = p_redf xr = p_xr xrud = p_xrud c c********************************************************************* c Mach number < 0.1 c********************************************************************* c if(mach.lt.0.1d0) goto 90 c c--- Additional constants c machq = mach*mach beta2 = 1.d0 - machq beta = sqrt(beta2) gamma = redf/beta2 mue = redf*mach/beta2 alpha = redf*mach*mach/beta2 c c--- Fixed node number for collacation method c n = 21 nint = 4*n + 1 c pin = pi/(2*n+1) pinint = pi/(2*nint+1) c ckon1 = cmplx(-0.25d0*beta,0.d0) ckon2 = (0.125d0*redf/beta)*ci c d2log = 1.d0 - 2.d0*log(2.d0) c c*********************************************************************** c Improper integral from minus infinity to minus 2 c for preparing the improper integral of the kernel function c*********************************************************************** c call unint(cu0) c c*********************************************************************** c 1. Computing the weights for integration of the non-singular part c of the kernel function c 2. Computing the coefficients of Jacobi's polynomials c for integration of the non-singular part of the kernel function c*********************************************************************** c do 5 m = 1,nint pm = 2*m*pinint gew(m) = 1.d0 - cos(pm) 5 continue c do 10 k = 1,n do 10 m = 1,nint pm = 2*m*pinint pmh = 0.5d0*pm fmjac(k,m) = sin((2.d0*k-1.d0)*pmh)/sin(pmh) 10 continue c c*********************************************************************** c Computing the kernel function at the abscissas (x,xsi) c*********************************************************************** c do 20 j = 1,n pj = (2*j-1)*pin pjh = 0.5d0*pj x = cos(pj) c c*********************************************************************** c Storing the required values of the kernel function c*********************************************************************** c do 25 m = 1,nint pm = 2*m*pinint pmh = 0.5d0*pm xsi = cos(pm) call kern(x,xsi,ckern) cvek(m) = ckern 25 continue c do 30 k = 1,n c dnenn = cos(pjh) fkjac = cos((2.d0*k-1.d0)*pjh)/dnenn fkjacm = cos((2.d0*k-3.d0)*pjh)/dnenn fkjacp = cos((2.d0*k+1.d0)*pjh)/dnenn c if(k.eq.1) then csing = ckon1 + ckon2*(fkjacp + d2log) else csing = ckon1*fkjac & + ckon2*((fkjacp+fkjac)/k - (fkjac+fkjacm)/(k-1)) end if c cint = (0.d0,0.d0) do 40 m = 1,nint cint = cint + gew(m)*cvek(m)*fmjac(k,m) 40 continue cint = 2.d0*pinint*cint c cmat(j,k) = csing + cint c 30 continue c 20 continue c 90 continue c open(16,file=outdat) c c################################################################ c Loop for computing pressure distributions and coefficients c for activated degrees of freedom c################################################################ c do 100 ifr = 1,3 c if(ifvek(ifr).ne.1) goto 100 c c c********************************************************************* c Computing coefficients for Mach numbers < 0.1 c********************************************************************* c if(mach.lt.0.1d0) then c call theo(ifr,redf,xr,xrud,cl,cmp,cmr) goto 120 c end if c c********************************************************************* c Prescribed wake function for plunge, pitch and aileron pitch c********************************************************************* c do 50 j = 1,n c if(ifr.eq.1) cw(j) = ci*redf c if(ifr.eq.2) then pj = (2*j-1)*pin x = cos(pj) cw(j) = -1.d0 - ci*redf*(x-xr) end if c if(ifr.eq.3) then pj = (2*j-1)*pin x = cos(pj) if(x.lt.xrud) then cw(j) = (0.d0,0.d0) else cw(j) = -1.d0 - ci*redf*(x-xrud) end if end if c 50 continue c c******************************************************** c Solution of the system of collocation equations c******************************************************** c call lingln(n,cmat,cw,ckof,cdet) c c********************************************************* c Computing and printing the pressure distribution c********************************************************* c if(idr.eq.1) then write(6,101) mach write(6,102) redf write(16,101) mach write(16,102) redf if(ifr.eq.1) write(6,103) if(ifr.eq.2) write(6,104) xr if(ifr.eq.3) write(6,105) xrud if(ifr.eq.1) write(16,103) if(ifr.eq.2) write(16,104) xr if(ifr.eq.3) write(16,105) xrud c 101 format(/10x,'Mach number',7x,' : ',f6.4) 102 format(/10x,'Reduced frequency : ',f6.4) 103 format(/10x,'DOF profile plunge') 104 format(/10x,'DOF profile pitch at xr = ',f6.3) 105 format(/10x,'DOF aileron pitch at xrud = ',f6.3) c write(6,61) write(16,61) 61 format(/10x,'Pressure distribution '/) c write(6,62) write(6,63) write(16,62) write(16,63) 62 format(14x,'x',9x,'Re(dcp)',7x,'Im(dcp)') 63 format(10x,38('-')) c do 60 j = 1,n pj = (2*j-1)*pin pjh = 0.5d0*pj x = cos(pj) wurz = sqrt((1.d0-x)/(1.d0+x)) c cdr = (0.d0,0.d0) do 65 k = 1,n rk = float(k) cdr = cdr + ckof(k)*sin((2.d0*rk-1.d0)*pjh)/sin(pjh) 65 continue cdr = wurz*cdr c write(6,64) x,cdr write(16,64) x,cdr 64 format(10x,f8.4,2(2x,f12.8)) c 60 continue c end if c c*************************************************** c Aerodynamic coefficients computed from c pressure distribution c*************************************************** c cl = pih*ckof(1) cmp = - 0.25d0*pi*(-ckof(1)*(0.5d0+xr) + 0.5d0*ckof(2)) c phic = acos(xrud) sinp = sin(phic) sin2p = sin(2.d0*phic) sin3p = sin(3.d0*phic) c cmr = ( - (0.5d0+xrud)*phic + (1.d0+xrud)*sinp & - 0.25d0*sin2p )*ckof(1) + ( 0.5d0*phic - (0.5d0+xrud)*sinp & + 0.5d0*(0.5d0+xrud)*sin2p - sin3p/6.d0 )*ckof(2) do 70 k = 3,n rk = float(k) rkm2 = rk - 2.d0 rkm1 = rk - 1.d0 rkp1 = rk + 1.d0 cmr = cmr + ( 0.5d0*sin(rkm2*phic)/rkm2 & - (0.5d0+xrud)*sin(rkm1*phic)/rkm1 & + (0.5d0+xrud)*sin(rk*phic)/rk & - 0.5d0*sin(rkp1*phic)/rkp1 )*ckof(k) 70 continue cmr = - 0.25d0*cmr c 120 continue c *** Compatibility with FLT2DQ: c 1. Plunge starts at the bottom (column 1) c 2. All moments mathematical negative oriented (rows 2,3) c Means for the signs c 1.- + + 2.- + + c - + + -> + - - c - + + + - - c for ifr = 1 2 3 c >>> Storing the data for common block POSSIO If (ifr .eq. 1) then CQVAL(1,ifr) = -CL CQVAL(2,ifr) = +CMP CQVAL(3,ifr) = +CMR else if (ifr .ge. 2) then CQVAL(1,ifr) = +CL CQVAL(2,ifr) = -CMP CQVAL(3,ifr) = -CMR end if 100 continue c c################################################################ c End of loop for pressure and aerodynamic coefficients c of active DOFs c################################################################ c close(16,status='keep') c c Print *,1,CQVAL(1,1),CQVAL(2,1),CQVAL(3,1) c Print *,2,CQVAL(1,2),CQVAL(2,2),CQVAL(3,2) c Print *,3,CQVAL(1,3),CQVAL(2,3),CQVAL(3,3) return end c subroutine kern(x,xsi,ckern) c***************************************************************** c Computing the non-singular part of the kernel function for c Possio's integral equation c***************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) real*8 mach,machq,mue c common /piwert/ pih,pi common /imag/ ci common /ma/ mach,machq,beta,beta2 common /rf/ redf common /gam/ mue,gamma,alpha common /uint/ cu0 c external cfunc3 c save c dif = x - xsi c write(6,11) dif 11 format(10x,'dif = ',f16.8) arg = mue*abs(dif) pifak = 2.d0/pi sig = dif/abs(dif) c c***************************************************************** c Terms which contain the Hankel functions H02 and H12 c***************************************************************** c call hank02(arg,ch0) call hank12(arg,ch1) c c***************************************************************** c Computing the remainder of the improper integral c***************************************************************** c if(dif.lt.0.d0) then call romint(-2.d0,dif,cfunc3,1.d-8,cu1) else call romint(-2.d0,-dif,cfunc3,1.d-8,cu1) call romint(-dif,dif,cfunc3,1.d-8,cu2) cu1 = cu1 + cu2 end if c ckern = exp(ci*alpha*dif)*(ci*machq*ch0 + mach*sig*ch1) & - exp(- ci*redf*dif)*( ci*redf*mach*(cu0+cu1) & - pifak*beta2*(cexint(2.d0*ci*gamma) & - cexint(-ci*gamma*dif) + ci*pih*(1.d0+sig)) ) c ckern = ckern - (ci*pifak/gamma)/dif - pifak*log(abs(dif)) c ckern = (ci*redf/(8.d0*beta))*ckern c return end c subroutine unint(cu0) c****************************************************************** c Computing the improper integral c c - 2 c - c | c | sign(t)*H12(mue*abs(t))*exp(i*gamma*t) dt c | c - c - infinity c c c parameter : c cu0 (output,komplex) : value of the integral c c****************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) real*8 mach,machq,mue c common /piwert/ pih,pi common /imag/ ci common /ma/ mach,machq,beta,beta2 common /rf/ redf common /gam/ mue,gamma,alpha common /euz/ euzahl c external cfunc1,cfunc2 c save c xtrenn = - 5000.d0 fak = redf/(1.d0-mach) xta = sqrt(fak*abs(xtrenn)) cfak1 = (2.d0/sqrt(pi))*(1.d0-ci)* & sqrt((1.d0-mach)*beta2/(redf*redf*mach)) c call romint(xtrenn,-2.d0,cfunc1,1.d-8,cinttr) call romint(0.d0,xta,cfunc2,1.d-8,cint0) cu0 = cfak1*(0.5d0*sqrt(0.5d0*pi)*(1.d0-ci) - cint0) + cinttr c return end c subroutine romint(a,b,cf,eps,cint) c****************************************************************** c Computing an integral of a complex function using c Romberg integration c Parameters: c a (input) : lower boundary c b (input) : upper boundary c cf : function which is integrated, in the calling c program declared as external c eps : given limit of accuracy c cint (output) : complex value of the integral c****************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) c dimension cmat(50,50) c c*************************************************** c Checking the integral's boundaries c*************************************************** c if(abs(b-a).lt.1.d-10) then cint = (0.d0,0.d0) goto 200 end if c c*************************************************** c Initial value c*************************************************** c intval = 32 h = (b-a)/intval j = 1 ke = intval/2 csum = (0.d0,0.d0) do 5 i = 1,intval+1 x = a + (i-1)*h csum = csum + cf(x) 5 continue cmat(j,1) = h*(csum - 0.5d0*cf(a) - 0.5d0*cf(b)) c write(6,31) intval,ke,cmat(j,1) c c***************************************************************** c Computing the integral with half the sampling points c using the function values already evaluated c***************************************************************** c 10 continue h = 0.5d0*h intval = 2*intval j = j+1 ke = 2*ke c csum = (0.d0,0.d0) do 20 i = 1,ke x = a + (2*i-1)*h csum = csum + cf(x) 20 continue cmat(j,1) = 0.5d0*cmat(j-1,1) + h*csum c c************************************************* c Romberg extrapolation c************************************************* c do 30 i = 2,j fak1 = 4**(i-1) fak2 = 1.d0/(fak1- 1.d0) cmat(j,i) = fak2*(fak1*cmat(j,i-1)-cmat(j-1,i-1)) 30 continue c c write(6,31) intval,ke,cmat(j,1),cmat(j,j) 31 format(2x,i8,2x,i8,4(2x,f14.8)) c c****************************************************************** c Stopping the computing after reaching the accuracy c****************************************************************** c if(j.gt.1) then c crel = (cmat(j,j)-cmat(j-1,j-1))/cmat(j-1,j-1) if(abs(crel).lt.eps) then goto 100 else goto 10 end if c else c goto 10 c end if c 100 continue cint = cmat(j,j) 200 return end c function cexint(cz) c****************************************************************** c Computing the complex exponential integral c c infinity c - c | exp(-t) c E(z) = |-------- dt c | t c - c z c c c parameters: c cz (input,komplex) : Argument for which computation is done c****************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) c common /euz/ euzahl c save c cexint = - euzahl - log(cz) c n = 1 can1 = - cz csum = can1 c 10 continue n = n + 1 can = - (n-1)*cz*can1/(n*n) csum = csum + can sum = abs(csum) c write(6,11) n,sum 11 format(10x,i3,2x,d12.4) can1 = can if(abs(can).lt.1.d-10) then goto 100 else goto 10 end if c 100 cexint = cexint - csum c return end c subroutine hank02(x,cf) c****************************************************************** c Computing the Hankel funktion 0. order and 2. kind H02(x) c parameters : c x (input,real) : sampling point for which h02 is computed c x must be grater than zero c cf (output,complex) : value H02 at x c****************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) c common /piwert/ pih,pi common /euz/ euzahl c save c if(x.le.0.d0) then write(6,1) x 1 format(2x,'x = ', d12.4,2x,'No computation of H02(x)!') goto 100 end if c eps = 1.d-12 pifak = 2.d0/pi ci = cmplx(0.d0,1.d0) c if(x.le.15.d0) then c ak = 1.d0 cbk = 1.d0 - ci*pifak*(log(0.5d0*x) + euzahl) csum = ak*cbk k = 0 c 10 continue k = k + 1 c fak = - 0.25d0*x*x/(k*k) ak = fak*ak cbk = cbk + ci*pifak/k csum = csum + ak*cbk c if(abs(ak*cbk).lt.eps) then cf = csum goto 100 else goto 10 end if c else c c---- Asymptotic approximation for large arguments c ak = 1.d0 sump = ak bk = - 1.d0/(8*x) sumq = bk k = 0 c 20 continue k = k + 1 c fakak = - (4*k-3)*(4*k-3)*(4*k-1)*(4*k-1)/((2*k-1)*2*k*8*x*8*x) ak = fakak*ak sump = sump + ak c fakbk = - (4*k-1)*(4*k-1)*(4*k+1)*(4*k+1)/((2*k+1)*2*k*8*x*8*x) bk = fakbk*bk sumq = sumq + bk c if(abs(ak).lt.eps.and.abs(bk).lt.eps) then arg = x - 0.25d0*pi cf = sqrt(pifak/x)*(sump - ci*sumq)*exp(-ci*arg) goto 100 else goto 20 end if c end if c 100 continue return end c subroutine hank12(x,cf) c****************************************************************** c Computing the Hankel funktion 1. order and 2. kind H02(x) c parameters : c x (input,real) : sampling point for which H12 is computed c x must be grater than zero c cf (output,complex) : value H12 at x c****************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) c common /piwert/ pih,pi common /euz/ euzahl c save c if(x.le.0.d0) then write(6,1) x 1 format(2x,'x = ', d12.4,2x,'No computation of H12(x)!') goto 100 end if c eps = 1.d-12 pifak = 2.d0/pi ci = cmplx(0.d0,1.d0) c if(x.le.15.d0) then c ak = 1.d0 cbk = 1.d0 - ci*pifak*(log(0.5d0*x) + euzahl - 0.5d0) csum = ak*cbk k = 0 c 10 continue k = k + 1 c fak = - 0.25d0*x*x/(k*(k+1)) ak = fak*ak cbk = cbk + ci*0.5d0*pifak*(1.d0/k + 1.d0/(k+1)) csum = csum + ak*cbk c if(abs(ak*cbk).lt.eps) then cf = ci*pifak/x + 0.5d0*x*csum goto 100 else goto 10 end if c else c c---- Asymptotic approximation for large arguments c ak = 1.d0 sump = ak bk = 3.d0/(8*x) sumq = bk k = 0 c 20 continue k = k + 1 c fakak = - (4-(4*k-3)*(4*k-3))*(4-(4*k-1)*(4*k-1)) & /((2*k-1)*2*k*8*x*8*x) ak = fakak*ak sump = sump + ak c fakbk = - (4-(4*k-1)*(4*k-1))*(4-(4*k+1)*(4*k+1)) & /((2*k+1)*2*k*8*x*8*x) bk = fakbk*bk sumq = sumq + bk c if(abs(ak).lt.eps.and.abs(bk).lt.eps) then arg = x - 0.75d0*pi cf = sqrt(pifak/x)*(sump - ci*sumq)*exp(-ci*arg) goto 100 else goto 20 end if c end if c 100 continue return end c function cfunc1(x) c**************************************************** c Computing the function c c F = sign(x)*H12(mue*abs(x))*exp(i*gamma*x) c c**************************************************** implicit real*8(a-b,d-h,o-z),complex*16(c) real*8 mue c common /piwert/ pih,pi common /imag/ ci common /gam/ mue,gamma,alpha c save c arg = mue*abs(x) sig = 1.d0 if(x.lt.0.d0) sig = - 1.d0 c call hank12(arg,ch) cfunc1 = sig*ch*exp(ci*gamma*x) c return end c function cfunc2(x) c**************************************************** c Computing the function c c F = exp(-i*x*x) c c**************************************************** implicit real*8(a-b,d-h,o-z),complex*16(c) c common /imag/ ci c save c cfunc2 = exp(-ci*x*x) c return end c function cfunc3(x) c*********************************************************************** c Computing the function c c F = sign(x)*( H12(mue*abs(x))-2.i/(pi*mue*abs(x)) )*exp(i*gamma*x) c c*********************************************************************** implicit real*8(a-b,d-h,o-z),complex*16(c) real*8 mue c common /piwert/ pih,pi common /imag/ ci common /gam/ mue,gamma,alpha c save c arg = mue*abs(x) sig = 1.d0 if(x.lt.0.d0) sig = - 1.d0 c if(abs(x).gt.1.d-12) then call hank12(arg,ch) cfunc3 = sig*(ch - 2.d0*ci/(pi*mue*abs(x)))*exp(ci*gamma*x) else cfunc3 = (0.d0,0.d0) end if c return end c subroutine lingln(n,ca,cb,cx,cdet) c*********************************************************************** c Computing a linear complex system of equations c ca*cx = cb with given n*n - matrix ca given c n-dimensional vector cb with the solution vector cx c*********************************************************************** c implicit real*8(a-b,d-h,o-z),complex*16(c) integer n,iz c parameter(nx=51,nxp1=nx+1) c dimension ca(nx,nx),cx(nx),cb(nx),caa(nx,nxp1),iz(nxp1) c n1=n+1 c do 1 i=1,n iz(i)=i caa(i,n1)=cb(i) do 1 k=1,n 1 caa(i,k)=ca(i,k) c cdet=(1.d0,0.d0) c do 5 j=1,n amax=0 c do 2 i=j,n aij=abs(caa(iz(i),j)) if(amax.lt.aij) imax=i 2 if(amax.lt.aij) amax=aij c jj=iz(imax) iz(imax)=iz(j) iz(j)=jj caaj=caa(jj,j) cdet=cdet*caaj c if(abs(caaj).eq.0.d0) go to 9 c j1=j+1 do 3 k=j1,n1 caa(jj,k)=caa(jj,k)/caaj if(abs(caa(jj,k)).lt.1.d-14) caa(jj,k) = (0.d0,0.d0) 3 continue c do 5 i=1,n if(i.eq.jj) go to 5 do 4 k=j1,n1 4 caa(i,k)=caa(i,k)-caa(i,j)*caa(jj,k) c 5 continue c do 6 i=1,n 6 cx(i)=caa(iz(i),n1) c if(n.eq.1) goto 9 do 7 j=2,n do 7 k=j,n 7 if(iz(j-1).gt.iz(k)) cdet=-cdet c 9 return end c subroutine theo(ifr,redf,xr,xrud,cl,cmp,cmr) c************************************************************************ c Analytic computation of aerodynamic coefficients for profile lift, c moment and aileron moment using Theordorsen's function c************************************************************************ c implicit real*8(a-b,d-h,o-z),complex*16(c) c common /piwert/ pih,pi common /imag/ ci c save c call hank02(redf,ch0) call hank12(redf,ch1) c ctheo = ch1/(ch1 + ci*ch0) c if(ifr.eq.1) then c cl = pi*redf*redf*( 1.d0 - ci*(2.d0/redf)*ctheo ) c cmp = pih*redf*redf*(xr - ci*(2.d0/redf)*(0.5d0+xr)*ctheo) c phic = acos(xrud) cmr = redf*redf*( 0.5d0*xrud*phic - 0.25d0*sin(phic) & - 0.25d0*xrud*sin(2.d0*phic) + sin(3.d0*phic)/12.d0 ) & - ci*redf*ctheo*( (0.5d0+xrud)*phic - (1.d0+xrud)*sin(phic) & + 0.25d0*sin(2.d0*phic) ) c end if c if(ifr.eq.2) then c phic = acos(xrud) sinp = sin(phic) sin2p = sin(2.d0*phic) sin3p = sin(3.d0*phic) sin4p = sin(4.d0*phic) c cl = pi*redf*redf*( xr + (2.d0/(redf*redf))*ctheo & + (ci/redf)*(1.d0 + (1.d0-2.d0*xr)*ctheo) ) c cmp = pih*redf*redf*(0.125d0 + xr*xr - ci/redf*(0.5d0-xr) & + (2.d0/redf)*(0.5d0+xr)*(1.d0/redf + ci*(0.5d0-xr))*ctheo) c cmr = redf*redf*( 0.5d0*phic*(0.125d0+xrud*xr) & - 0.25d0*(xr+0.5d0*xrud)*sinp & - 0.25d0*xrud*xr*sin2p & + (xr+0.5d0*xrud)*sin3p/12.d0 & - sin4p/64.d0 ) c & - ctheo*( - (0.5d0+xrud)*phic + (1.d0+xrud)*sinp & - 0.25d0*sin2p ) c & - ci*redf* ( 1.5d0*( -0.5d0*xrud*phic + 0.25d0*sinp & + 0.25d0*xrud*sin2p - sin3p/12.d0 ) c & + 0.5d0*( 0.5d0*(1.d0+xrud)*phic - (0.75d0+xrud)*sinp & + 0.25d0*(1.d0+xrud)*sin2p - sin3p/12.d0 ) c & + (0.5d0-xr)*ctheo*( - (0.5d0+xrud)*phic & + (1.d0+xrud)*sinp - 0.25d0*sin2p ) ) c end if c if(ifr.eq.3) then c phic = acos(xrud) sinp = sin(phic) phicq = phic*phic sinpq = sinp*sinp xrudq = xrud*xrud c dk1 = xrud*phic - (2.d0+xrudq)*sinp/3.d0 dk2 = 2.d0*(sinp+phic) dk3 = - (xrud*sinp-phic) dk4 = (1.d0-2.d0*xrud)*phic + (2.d0-xrud)*sinp c cl = redf*redf*dk1 + dk2*ctheo + ci*redf*(dk3 + dk4*ctheo) c dk1 = - 0.5d0*(xrud+1.d0)*sinp dk2 = (sinp+phic)*(xr+0.5d0) dk3 = (0.125d0 + xrudq)*phic & - 0.125d0*xrud*sinp*(7.d0 + 2.d0*xrudq) & + (xrud-xr)*(sinp*(2.d0+xrudq)/3.d0 - xrud*phic) dk4 = sinp*(xrudq-1.d0)/3.d0 - (xrud-xr)*(xrud*sinp-phic) dk5 = (1.d0-2.d0*xrud)*phic + (2.d0-xrud)*sinp c cmp = dk1 + dk2*ctheo + 0.5d0*redf*redf*dk3 & - 0.5d0*ci*redf*( dk4 + dk5*(0.5d0 - (xr+0.5d0)*ctheo)) c dk1 = (1.d0+xrud)*sinp*(phic-sinp) dk2 = 0.25d0*xrud*sinp*phic*(7.d0+2.d0*xrudq) & - (0.125d0+xrudq)*phicq & - 0.125d0*(1.d0-xrudq)*(5.d0*xrudq+4.d0) dk3 = (2.d0+xrud)*sinpq + (1.d0-xrud)*phic*sinp & - (2.d0*xrud+1.d0)*phicq c dk4 = (0.5d0-xrud)*phicq + (xrudq-xrud+1.d0)*phic*sinp & - (1.d0-0.5d0*xrud)*xrud*sinpq dk5 = 2.d0*(xrudq-0.25d0)*phicq - 3.d0*xrud*phic*sinp & + 2.d0*(1.d0-0.25d0*xrudq)*sinpq c cmr = dk1 + redf*redf*dk2 + dk3*ctheo & + ci*redf*( dk4 + dk5*ctheo ) cmr = - cmr/(2.d0*pi) c end if c return end