********************* Delete above including this line ***************** c------------------------------------------------------------------ c file = "twodd_new.f". Maximum allowable boundary elements = 200 c Dimension statement changed by Jian Lin on Nov. 22, 1998 c c------------------------------------------------------------------ c program twodd: 2-d displacement discontinuity method c c------------------------------------------------------------------ character*16 ifile,ofile dimension xm(200),ym(200),a(200),cosbet(200),sinbet(200), 1 lmntype(200),kod(200) c complex cmplx,conjg,csqrt,clog,cexp complex gams,gamn,gamsbar,gamnbar,ckappa common /s1/ pi,pr,pr1,pr2,con,cons,rkappa common /s2/ sxxs,sxxn,syys,syyn,sxys,sxyn,uxs,uxn,uys,uyn common /s3/ c(400,400),b(400),d(400) common /s4/ gams,gamn,gamsbar,gamnbar,ckappa c c open input/output files c write(*,*) 'input file is?' read(*,*) ifile write(*,*) 'output file is?' read(*,*) ofile open (unit=5,file=ifile) open (unit=6,file=ofile) c read(5,*) numbs,numos,ksym read(5,*) pr,e,xsym,ysym read(5,*) pxx,pyy,pxy write(6,6) numbs,numos go to (80,85,90,95),ksym 80 write(6,7) go to 100 85 write(6,8) xsym go to 100 90 write(6,9) ysym go to 100 95 write(6,10) xsym,ysym c 100 continue write(6,11) pr,e write(6,12) pxx,pyy,pxy c pi = 4.*atan(1.) con = 1./(4.*pi*(1.-pr)) cons = e/(1.+pr) pr1 = 1.-2.*pr pr2 = 2.*(1.-pr) c c gams = constant for s unit displacement c gamn = constant for n unit displacement c rkappa = 3.-4.*pr gamma = -cons*con/2. gams = cmplx(0.,-gamma) gamn = cmplx(gamma,0.) gamsbar = conjg(gams) gamnbar = conjg(gamn) ckappa = cmplx(rkappa,0.) c c define locations, sizes, orientations and boundary conditions of c boundary elements c numbe = 0 do 110 n=1,numbs read(5,*) num,xbeg,ybeg,xend,yend,kode,bvs,bvn,ltype xd = (xend-xbeg)/num yd = (yend-ybeg)/num sw = sqrt(xd*xd+yd*yd) c do 110 ne = 1,num numbe = numbe + 1 m = numbe xm(m) = xbeg + 0.5*(2.*ne-1.)*xd ym(m) = ybeg + 0.5*(2.*ne-1.)*yd a(m) = 0.5*sw sinbet(m) = yd/sw cosbet(m) = xd/sw kod(m) = kode lmntype(m) = ltype mn = 2*m ms = mn - 1 b(ms) = bvs b(mn) = bvn 110 continue c write(6,13) do 115 m=1,numbe size = 2.*a(m) angle = 180.*atan2(sinbet(m),cosbet(m))/pi write(6,15) m,kod(m),xm(m),ym(m),size,angle, 1 b(2*m-1),b(2*m),lmntype(m) 115 continue c c adjust stress boundary values to account for (uniform) initial stresses c do 150 n = 1,numbe nn = 2*n ns = nn - 1 cosb = cosbet(n) sinb = sinbet(n) sigs = (pyy-pxx)*sinb*cosb + pxy*(cosb*cosb-sinb*sinb) sign = pxx*sinb*sinb - 2.*pxy*sinb*cosb + pyy*cosb*cosb go to (120,150,130,140),kod(n) 120 b(ns) = b(ns) - sigs b(nn) = b(nn) - sign go to 150 130 b(nn) = b(nn) - sign go to 150 140 b(ns) = b(ns) - sigs 150 continue c c compute influence coefficients and set up system of c algebraic equations c do 300 i = 1,numbe in = 2*i is = in - 1 xi = xm(i) yi = ym(i) cosbi = cosbet(i) sinbi = sinbet(i) kode = kod(i) c do 300 j = 1,numbe jn = 2*j js = jn - 1 call initl xj = xm(j) yj = ym(j) cosbj = cosbet(j) sinbj = sinbet(j) aj = a(j) ltypj = lmntype(j) c call coeff(xi,yi,xj,yj,aj,cosbj,sinbj,+1,ltypj) c go to (240,210,220,230),ksym c 210 xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,cosbj,-sinbj,-1,ltypj) go to 240 c 220 yj = 2.*ysym - ym(j) call coeff(xi,yi,xj,yj,aj,-cosbj,sinbj,-1,ltypj) go to 240 c 230 xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,cosbj,-sinbj,-1,ltypj) xj = xm(j) yj = 2.*ysym - ym(j) call coeff(xi,yi,xj,yj,aj,-cosbj,sinbj,-1,ltypj) xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,-cosbj,-sinbj,+1,ltypj) c 240 continue c go to (250,260,270,280),kode c 250 c(is,js)= (syys-sxxs)*sinbi*cosbi+sxys*(cosbi*cosbi-sinbi*sinbi) c(is,jn)= (syyn-sxxn)*sinbi*cosbi+sxyn*(cosbi*cosbi-sinbi*sinbi) c(in,js)= sxxs*sinbi*sinbi-2.*sxys*sinbi*cosbi+syys*cosbi*cosbi c(in,jn)= sxxn*sinbi*sinbi-2.*sxyn*sinbi*cosbi+syyn*cosbi*cosbi go to 300 c 260 c(is,js)= uxs*cosbi+uys*sinbi c(is,jn)= uxn*cosbi+uyn*sinbi c(in,js)= -uxs*sinbi+uys*cosbi c(in,jn)= -uxn*sinbi+uyn*cosbi go to 300 c 270 c(is,js)= uxs*cosbi+uys*sinbi c(is,jn)= uxn*cosbi+uyn*sinbi c(in,js)= sxxs*sinbi*sinbi-2.*sxys*sinbi*cosbi+syys*cosbi*cosbi c(in,jn)= sxxn*sinbi*sinbi-2.*sxyn*sinbi*cosbi+syyn*cosbi*cosbi go to 300 c 280 c(is,js)= (syys-sxxs)*sinbi*cosbi+sxys*(cosbi*cosbi-sinbi*sinbi) c(is,jn)= (syyn-sxxn)*sinbi*cosbi+sxyn*(cosbi*cosbi-sinbi*sinbi) c(in,js)= -uxs*sinbi+uys*cosbi c(in,jn)= -uxn*sinbi+uyn*cosbi c 300 continue c c solve system of algebraic equations c n = 2*numbe call solve(n) c c compute boundary displacements and stresses c write(6,16) do 600 i=1,numbe in = 2*i is = in - 1 xi = xm(i) yi = ym(i) cosbi = cosbet(i) sinbi = sinbet(i) ltypi = lmntype(i) c uxneg = 0. uyneg = 0. sigxx = pxx sigyy = pyy sigxy = pxy c do 570 j = 1,numbe jn = 2*j js = jn - 1 call initl xj = xm(j) yj = ym(j) aj = a(j) cosbj = cosbet(j) sinbj = sinbet(j) ltypj = lmntype(j) c call coeff(xi,yi,xj,yj,aj,cosbj,sinbj,+1,ltypj) c go to (540,510,520,530),ksym c 510 xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,cosbj,-sinbj,-1,ltypj) go to 540 c 520 yj = 2.*ysym - ym(j) call coeff(xi,yi,xj,yj,aj,-cosbj,sinbj,-1,ltypj) go to 540 c 530 xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,cosbj,-sinbj,-1,ltypj) xj = xm(j) yj = 2.*ysym - ym(j) call coeff(xi,yi,xj,yj,aj,-cosbj,sinbj,-1,ltypj) xj = 2.*xsym - xm(j) call coeff(xi,yi,xj,yj,aj,-cosbj,-sinbj,+1,ltypj) c 540 continue c uxneg = uxneg + uxs*d(js) + uxn*d(jn) uyneg = uyneg + uys*d(js) + uyn*d(jn) sigxx = sigxx + sxxs*d(js) + sxxn*d(jn) sigyy = sigyy + syys*d(js) + syyn*d(jn) sigxy = sigxy + sxys*d(js) + sxyn*d(jn) c 570 continue c if(ltypi.eq.2) then uxneg = -uxneg uyneg = -uyneg endif c usneg = uxneg*cosbi + uyneg*sinbi unneg = -uxneg*sinbi + uyneg*cosbi uspos = usneg - d(is) unpos = unneg - d(in) uxpos = uspos*cosbi - unpos*sinbi uypos = uspos*sinbi + unpos*cosbi sigs = (sigyy-sigxx)*sinbi*cosbi+sigxy*(cosbi*cosbi-sinbi*sinbi) sign = sigxx*sinbi*sinbi-2.*sigxy*sinbi*cosbi+sigyy*cosbi*cosbi c write(6,17) i,d(is),usneg,uspos,d(in),unneg,unpos,uxneg,uyneg, 1 uxpos,uypos,sigs,sign c 600 continue c c compute displacements and stresses at specified points in body c if (numos.le.0) go to 900 write(6,18) npoint = 0 do 890 n=1,numos read(5,*) xbeg,ybeg,xend,yend,numpb nump = numpb + 1 delx = (xend-xbeg)/nump dely = (yend-ybeg)/nump if (numpb.gt.0) nump = nump+1 if (delx**2+dely**2.eq.0.) nump = 1 c do 890 ni = 1,nump xp = xbeg + (ni-1)*delx yp = ybeg + (ni-1)*dely c ux = 0. uy = 0. sigxx = pxx sigyy = pyy sigxy = pxy c do 880 j = 1,numbe jn = 2*j js = jn-1 call initl xj = xm(j) yj = ym(j) aj = a(j) c if (sqrt((xp-xj)**2+(yp-yj)**2).lt.2.*aj) go to 890 c cosbj = cosbet(j) sinbj = sinbet(j) ltypj = lmntype(j) c call coeff(xp,yp,xj,yj,aj,cosbj,sinbj,+1,ltypj) c go to (840,810,820,830),ksym c 810 xj = 2.*xsym - xm(j) call coeff(xp,yp,xj,yj,aj,cosbj,-sinbj,-1,ltypj) go to 840 c 820 yj = 2.*ysym - ym(j) call coeff(xp,yp,xj,yj,aj,-cosbj,sinbj,-1,ltypj) go to 840 c 830 xj = 2.*xsym - xm(j) call coeff(xp,yp,xj,yj,aj,cosbj,-sinbj,-1,ltypj) xj = xm(j) yj = 2.*ysym - ym(j) call coeff(xp,yp,xj,yj,aj,-cosbj,sinbj,-1,ltypj) xj = 2.*xsym - xm(j) call coeff(xp,yp,xj,yj,aj,-cosbj,-sinbj,+1,ltypj) c 840 continue c ux = ux + uxs*d(js) + uxn*d(jn) uy = uy + uys*d(js) + uyn*d(jn) sigxx = sigxx + sxxs*d(js) + sxxn*d(jn) sigyy = sigyy + syys*d(js) + syyn*d(jn) sigxy = sigxy + sxys*d(js) + sxyn*d(jn) c 880 continue c npoint = npoint + 1 write(6,20) npoint,xp,yp,ux,uy,sigxx,sigyy,sigxy c 890 continue c 900 continue c c find crack tip stress intensities (if crack tips are present) c call k_factor(lmntype,a,xm,ym,numbe) c c format statements c 2 format (/a80/) 6 format (/,109h number of straight-line segments (each containing a +t least one boundary element) used to define boundaries =,i3,//,12 +3h number of straight-line segments used to specify other location +s (i.e., not on a boundary) where results are to be found =,i3) 7 format (/,32h no symmetry conditions imposed.) 8 format (/,18h the line x = xs =,f12.4,23h is a line of symmetry.) 9 format (/,18h the line y = ys =,f12.4,23h is a line of symmetry.) 10 format (/,19h the lines x = xs =,f12.4,13h and y = ys =,f12.4, + 23h are lines of symmetry.) 11 format (/,18h poisson's ratio =,f6.2,//,18h young's modulus =, + e11.4) 12 format (/,31h xx-component of field stress =,e11.4,//,31h yy-compo +nent of field stress =,e11.4,//,31h xy-component of field stress = +,e11.4) 13 format (//,' boundary element data.',//, 1 ' element kode x (midpt) y (midpt) length', 2 ' angle us or sigma-s un or sigma-n lmntype',/) 15 format (2i6,1p3e12.4,1pe12.4,1p2e15.4,i8) 16 format (//,' displacements and stresses at midpoints', 1' of boundary elements.',//,' element ds us(-)', 2' us(+) dn un(-) un(+) ux(-) uy(-)', 3' ux(+) uy(+) sigma-s sigma-n',/) 17 format (i5,3x,1p12e10.2) 18 format (//,' displacements and stresses at specified', 1' points in the body.',//,' point x co-ord y co-ord', 2' ux uy sigxx sigyy sigxy',//) 20 format (i9,1p2e12.4,1p5e12.3) c stop end c c--------------------- subroutine initl c--------------------- common /s2/sxxs,sxxn,syys,syyn,sxys,sxyn,uxs,uxn,uys,uyn c sxxs = 0. sxxn = 0. syys = 0. syyn = 0. sxys = 0. sxyn = 0. c uxs = 0. uxn = 0. uys = 0. uyn = 0. c return end c c------------------------ subroutine solve(n) c------------------------ common /s3/ a(400,400),b(400),x(400) c nb = n - 1 do 20 j=1,nb l = j + 1 do 20 jj=l,n xm = a(jj,j)/a(j,j) do 10 i=j,n 10 a(jj,i) = a(jj,i) - a(j,i)*xm 20 b(jj) = b(jj) - b(j)*xm c x(n) = b(n)/a(n,n) do 40 j=1,nb jj = n-j l = jj + 1 sum = 0. do 30 i=l,n 30 sum = sum + a(jj,i)*x(i) 40 x(jj) = (b(jj)-sum)/a(jj,jj) c return end c c----------------------------------------------- subroutine k_factor(lmntype,a,xm,ym,numbe) c----------------------------------------------- c logical first c dimension lmntype(200),xm(200),ym(200),a(200) c common /s1/ pi,pr,pr1,pr2,con,cons,rkappa common /s3/ c(400,400),b(400),d(400) c first = .true. do 100 i=1,numbe in = 2*i is = in - 1 xi = xm(i) yi = ym(i) ai = a(i) ltypi = lmntype(i) c if(ltypi.eq.1 .or.ltypi.eq.2) then c if(first) then write(6,1) 1 format(///,' crack tip stress intensity factor(s)',/) first = .false. endif c sif1 = -cons*con*pi*sqrt(pi/2.)*d(in)/sqrt(ai) sif2 = -cons*con*pi*sqrt(pi/2.)*d(is)/sqrt(ai) theta = growdir(sif1,sif2)*180./pi write(6,2) xi,yi,i,sif1,sif2,theta 2 format(' for crack tip centered at x,y =', 1 1pe10.3,',',1pe10.3,' (element',i3,') ki,kii are ',1pe10.3, 2 ',',1pe10.3,' growth angle = ',0pf7.2/) endif c 100 continue c return end c c-------------------------------- function growdir(sif1,sif2) c-------------------------------- c c solve for takeoff angle of crack extension c by newton-rapheson solution of c k1*sin(growdir) + k2*(3*cos(growdir)-1) = 0 c thet = 0.0 s21 = sif2/sif1 10 continue sinthet = sin(thet) costhet = cos(thet) h = (s21*(3.*costhet-1.) + sinthet)/ 1 (3.*s21*sinthet - costhet) thet = thet + h if(abs(h).gt.1.0e-4) go to 10 c c..if not well converged then reiterate c growdir = thet c return end c c------------------------------------------------------- subroutine coeff(x,y,cx,cy,a,cosb,sinb,msym,ltype) c------------------------------------------------------- if(ltype.eq.0) then c..constant displacement element call coeff0(x,y,cx,cy,a,cosb,sinb,msym) else if(ltype.eq.1) then c..crack tip element -- tip at beginning of element call coeff1(x,y,cx,cy,a,cosb,sinb,msym) else if(ltype.eq.2) then c..crack tip element -- tip at end of element call coeff2(x,y,cx,cy,a,cosb,sinb,msym) endif c return end c c-------------------------------------------------- subroutine coeff0(x,y,cx,cy,a,cosb,sinb,msym) c-------------------------------------------------- complex cmplx,conjg,csqrt,clog,cexp complex zb,zbbar,gams,gamn,gamsbar,gamnbar,expib,expmi2b, 1 zma,zpa,zma2,zpa2,zbarma,zbarpa,rzma,rzpa,ckappa,ca, 2 capphi,phi,psi,phiprime,phi2prime,psiprime, 3 cusn,cuxy,clnz c common /s1/ pi,pr,pr1,pr2,con,cons,rkappa common /s2/ sxxs,sxxn,syys,syyn,sxys,sxyn,uxs,uxn,uys,uyn common /s4/ gams,gamn,gamsbar,gamnbar,ckappa c cos2b = cosb*cosb - sinb*sinb sin2b = 2.*sinb*cosb expib = cmplx(cosb,sinb) expmi2b = cmplx(cos2b,-sin2b) cosb2 = cosb*cosb sinb2 = sinb*sinb c xb = (x-cx)*cosb + (y-cy)*sinb yb = -(x-cx)*sinb + (y-cy)*cosb zb = cmplx(xb,yb) if(zb.eq.(0.,0.)) zb = (0.,-0.000001) zbbar = conjg(zb) ca = cmplx(a,0.) zma = zb - ca zbarma = zbbar - ca zma2 = zma*zma zpa = zb + ca zbarpa = zbbar + ca zpa2 = zpa*zpa rzma = 1./zma rzpa = 1./zpa clnz = clog(zma) - clog(zpa) c phi = gams*clnz psi = gamsbar*clnz - gams*ca*(rzma+rzpa) phiprime = gams*(rzma-rzpa) capthet = 2.*real(gams*(rzma-rzpa)) capphi = gamsbar*(rzma-rzpa) - gams*(zbarma/zma2-zbarpa/zpa2) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxds = real(cuxy)/cons uyds = aimag(cuxy)/cons c sxxds = capthet - real(capphi) syyds = capthet + real(capphi) sxyds = aimag(capphi) phi = gamn*clnz psi = gamnbar*clnz - gamn*ca*(rzma+rzpa) phiprime = gamn*(rzma-rzpa) capthet = 2.*real(gamn*(rzma-rzpa)) capphi = gamnbar*(rzma-rzpa) - gamn*(zbarma/zma2-zbarpa/zpa2) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxdn = real(cuxy)/cons uydn = aimag(cuxy)/cons c sxxdn = capthet - real(capphi) syydn = capthet + real(capphi) sxydn = aimag(capphi) c uxs = uxs + msym*uxds uxn = uxn + uxdn uys = uys + msym*uyds uyn = uyn + uydn c sxxs = sxxs + msym*sxxds sxxn = sxxn + sxxdn syys = syys + msym*syyds syyn = syyn + syydn sxys = sxys + msym*sxyds sxyn = sxyn + sxydn c return end c c-------------------------------------------------- subroutine coeff1(x,y,cx,cy,a,cosb,sinb,msym) c-------------------------------------------------- complex cmplx,conjg,csqrt,clog,cexp complex zb,zbbar,gams,gamn,gamsbar,gamnbar,expib,expmi2b, 1 ckappa,c2a,clnz,capphi,phi,psi,phiprime,cusn,cuxy c common /s1/ pi,pr,pr1,pr2,con,cons,rkappa common /s2/ sxxs,sxxn,syys,syyn,sxys,sxyn,uxs,uxn,uys,uyn common /s4/ gams,gamn,gamsbar,gamnbar,ckappa c cos2b = cosb*cosb - sinb*sinb sin2b = 2.*sinb*cosb expib = cmplx(cosb,sinb) expmi2b = cmplx(cos2b,-sin2b) cosb2 = cosb*cosb sinb2 = sinb*sinb c xb = (x-cx)*cosb + (y-cy)*sinb + a yb = -(x-cx)*sinb + (y-cy)*cosb zb = cmplx(xb,yb) zbbar = conjg(zb) c2a = cmplx(2.*a,0.) clnz = clog((csqrt(zb)+csqrt(c2a))/(csqrt(zb)-csqrt(c2a))) sqr2 = sqrt(2.) c phi = -gams*csqrt(zb/a)*clnz phiprime = -gams/(csqrt(4.*a*zb))*clnz + sqr2*gams/(zb-2.*a) psi = (gams/2.-gamsbar)*csqrt(zb/a)*clnz 1 - 2.*sqr2*gams*a/(zb-2.*a) c capthet = 2.*real(phiprime) capphi = (gams+gams*zbbar/zb-2.*gamsbar)/(4.*csqrt(a*zb))*clnz 1 - (gams/sqr2*(1.-zbbar/zb) 2 + sqr2*(gams*(zbbar-2.*a)/(zb-2.*a) 3 - gamsbar))/(zb-2.*a) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxds = real(cuxy)/cons uyds = aimag(cuxy)/cons c sxxds = capthet - real(capphi) syyds = capthet + real(capphi) sxyds = aimag(capphi) c phi = -gamn*csqrt(zb/a)*clnz phiprime = -gamn/(csqrt(4.*a*zb))*clnz + sqr2*gamn/(zb-2.*a) psi = (gamn/2.-gamnbar)*csqrt(zb/a)*clnz 1 - 2.*sqr2*gamn*a/(zb-2.*a) c capthet = 2.*real(phiprime) capphi = (gamn+gamn*zbbar/zb-2.*gamnbar)/(4.*csqrt(a*zb))*clnz 1 - (gamn/sqr2*(1.-zbbar/zb) 2 + sqr2*(gamn*(zbbar-2.*a)/(zb-2.*a) 3 - gamnbar))/(zb-2.*a) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxdn = real(cuxy)/cons uydn = aimag(cuxy)/cons c sxxdn = capthet - real(capphi) syydn = capthet + real(capphi) sxydn = aimag(capphi) c uxs = uxs + msym*uxds uxn = uxn + uxdn uys = uys + msym*uyds uyn = uyn + uydn c sxxs = sxxs + msym*sxxds sxxn = sxxn + sxxdn syys = syys + msym*syyds syyn = syyn + syydn sxys = sxys + msym*sxyds sxyn = sxyn + sxydn c return end c c-------------------------------------------------- subroutine coeff2(x,y,cx,cy,a,cosb,sinb,msym) c-------------------------------------------------- complex cmplx,conjg,csqrt,clog,cexp complex zb,zbbar,gams,gamn,gamsbar,gamnbar,expib,expmi2b, 1 ckappa,c2a,clnz,capphi,phi,psi,phiprime,cusn,cuxy c common /s1/ pi,pr,pr1,pr2,con,cons,rkappa common /s2/ sxxs,sxxn,syys,syyn,sxys,sxyn,uxs,uxn,uys,uyn common /s4/ gams,gamn,gamsbar,gamnbar,ckappa c cos2b = cosb*cosb - sinb*sinb sin2b = 2.*sinb*cosb expib = cmplx(cosb,sinb) expmi2b = cmplx(cos2b,-sin2b) cosb2 = cosb*cosb sinb2 = sinb*sinb c xb = (x-cx)*cosb + (y-cy)*sinb - a yb = -(x-cx)*sinb + (y-cy)*cosb zb = cmplx(xb,yb) zbbar = conjg(zb) c2a = cmplx(2.*a,0.) clnz = clog((csqrt(-zb)+csqrt(c2a))/(csqrt(-zb)-csqrt(c2a))) sqr2 = sqrt(2.) c phi = gams*csqrt(-zb/a)*clnz phiprime = -gams/(csqrt(-4.*a*zb))*clnz - sqr2*gams/(zb+2.*a) c psi = (gamsbar-gams/2.)*csqrt(-zb/a)*clnz 1 - 2.*sqr2*gams*a/(zb+2.*a) c capthet = 2.*real(phiprime) capphi = (gams+gams*zbbar/zb-2.*gamsbar)/(4.*csqrt(-a*zb))*clnz 1 + (gams/sqr2*(1.-zbbar/zb) 2 + sqr2*gams*(zbbar+2.*a)/(zb+2.*a) 3 - sqr2*gamsbar)/(zb+2.*a) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxds = real(cuxy)/cons uyds = aimag(cuxy)/cons c sxxds = capthet - real(capphi) syyds = capthet + real(capphi) sxyds = aimag(capphi) c phi = gamn*csqrt(-zb/a)*clnz phiprime = -gamn/(csqrt(-4.*a*zb))*clnz - sqr2*gamn/(zb+2.*a) c psi = (gamnbar-gamn/2.)*csqrt(-zb/a)*clnz 1 - 2.*sqr2*gamn*a/(zb+2.*a) c capthet = 2.*real(phiprime) capphi = (gamn+gamn*zbbar/zb-2.*gamnbar)/(4.*csqrt(-a*zb))*clnz 1 + (gamn/sqr2*(1.-zbbar/zb) 2 + sqr2*gamn*(zbbar+2.*a)/(zb+2.*a) 3 - sqr2*gamnbar)/(zb+2.*a) c cusn = ckappa*phi - zb*conjg(phiprime) - conjg(psi) cuxy = cusn*expib c capphi = capphi*expmi2b c uxdn = real(cuxy)/cons uydn = aimag(cuxy)/cons c sxxdn = capthet - real(capphi) syydn = capthet + real(capphi) sxydn = aimag(capphi) c uxs = uxs + msym*uxds uxn = uxn + uxdn uys = uys + msym*uyds uyn = uyn + uydn c sxxs = sxxs + msym*sxxds sxxn = sxxn + sxxdn syys = syys + msym*syyds syyn = syyn + syydn sxys = sxys + msym*sxyds sxyn = sxyn + sxydn c c return end