***************************************************************** * qs2g20.f NATORI Hiroshi Yamanashi Univ. Dec. 1995 * ***************************************************************** implicit double complex(c) IMPLICIT double precision (A-B,D-H,O-Z) parameter(lx = 512) parameter(ly = 512) parameter(lxg = 512) parameter(lyg = 512) parameter(ione = 1) parameter(itwo = 2) parameter(ithr = 3) parameter(one = 1.0d0) parameter(two = 2.0d0) parameter(thr = 3.0d0) parameter(fr = 4.0d0) parameter(ten = 10.0d0) parameter(zero = 0.0d0) parameter(pi = 3.14159265d0) parameter(hlf = 0.5d0) parameter(t8 = 256.0d0) dimension cpsi(lx+5,ly+5), cev0(lx+5,ly+5) cone = cmplx(1.0d0,0.0d0) ci = cmplx(0.0d0,1.0d0) rone = 1.0d0 rzero = 0.0d0 rnzero = 0.02d0 tpi = 2.0d0*pi rt6 = sqrt(6.0d0) rt2 = sqrt(2.0d0) lxh = int(dble(lx)/2.0d0) lyh = int(dble(ly)/2.0d0) xmn = -10.0d0 xmx = 10.0d0 ymn = -10.0d0 ymx = 10.0d0 c ymn = -6.25d0 c ymx = 6.25d0 c eps = 0.5d0 eps = 0.09375d0 px = 3.14159265d0 py = 0.0d0 sigmax = 1.0d0 sigmay = 1.0d0 potr = 1.5d0 xa = -5.0d0 ya = 0.0d0 v = -pi*pi e = 1.0d0 n = 32500 k1 = 250 n = n+1 c ig = 6 c cieps = cmplx(0.0d0,-eps) cal = (cone+exp(cieps))*hlf cbt = (cone-exp(cieps))*hlf cgm = exp(cieps)*hlf dx = (xmx-xmn)/(lx) dy = (ymx-ymn)/(ly) dt = dx*dx*eps r1 = 8.0d0/(2.0d0*3.1415927d0* & sqrt(sigmax*sigmax*sigmay*sigmay)) rx2 = 2.0d0*sigmax*sigmax ry2 = 2.0d0*sigmay*sigmay m0 = 1 p = 1.0d0 **************************************************************** y = ymn do 23000 j = 3,ly+2 x = xmn do 23002 i = 3,lx+2 rtmp1 = -(x-xa)*(x-xa)/rx2-(y-ya)*(y-ya)/ry2 rtmp2 = px*x+py*y ctmp1 = cmplx(rtmp1,rtmp2) if ((mod(i,20).eq.1).and.(mod(j,20).eq.1)) then endif cpsi(i,j) = r1*exp(ctmp1) x = x+dx 23002 continue y = y+dy 23000 continue do 450 j=3,ly+2 cpsi(1,j)=cpsi(lx+1,j) cpsi(2,j)=cpsi(lx+2,j) cpsi(lx+3,j)=cpsi(3,j) cpsi(lx+4,j)=cpsi(4,j) 450 continue do 550 i=3,lx+2 cpsi(i,1)=cpsi(i,ly+1) cpsi(i,2)=cpsi(i,ly+2) cpsi(i,ly+3)=cpsi(i,3) cpsi(i,ly+4)=cpsi(i,4) 550 continue *************************************************************** y = ymn do 25570 j = 3,ly+2 x = xmn do 25575 i = 3,lx+2 if (x*x+y*y.le.potr*potr) then rtmp2 = v else rtmp2 = 0.0d0 endif c rtmp1 = (rtmp2+ten*sin(fr*pi*x)*sin(fr*pi*y)) rtmp1 = rtmp2 cev0(i,j) = exp( & cmplx(0.0d0,-rtmp1*dt) & ) x = x+dx 25575 continue y = y+dy 25570 continue do 60040 j = 3,ly+2 cev0(1,j) = cev0(lx+1,j) cev0(2,j) = cev0(lx+2,j) cev0(lx+3,j) = cev0(3,j) cev0(lx+4,j) = cev0(4,j) 60040 continue do 60050 i = 3,lx+2 cev0(i,1) = cev0(i,ly+1) cev0(i,2) = cev0(i,ly+2) cev0(i,ly+3) = cev0(i,3) cev0(i,ly+4) = cev0(i,4) 60050 continue *************************************************************** nfxv = 999 nf = 1999 do 53016 ka = 1,n if (mod(ka,k1).eq.1) then nf = nf+1 nfxv = nfxv+1 *************************************************************** open(unit = nfxv,status = "new") s = 0.0d0 write(nfxv,"(a)") "P3" write(nfxv,"(2i4)") lxg,lyg write(nfxv,"(a)") "255" write(nf,"(a)") "P2" write(nf,"(2i4)") lxg,lyg write(nf,"(a)") "255" do 23020 j = 3,ly+4 do 23022 i = 3,lx+2 psix = dble(cpsi(i,j)) psiy = imag(cpsi(i,j)) ra2 = psix*psix+psiy*psiy ira2 = int(512.0d0*ra2) if (ira2.gt.255) then ira2 = 255 endif s = s+ra2*dx*dy ra = sqrt(ra2) rtmp1 = hlf+sign(rone,ra2-rone) & *(hlf-ra/(rone+ra2)) rtmp2 = psix/(rt6*(rone+ra2)) rtmp3 = psiy/(rt2*(rone+ra2)) ir = int((rtmp1+two*rtmp2)*t8) ig = int((rtmp1-rtmp2+rtmp3)*t8) ib = int((rtmp1-rtmp2-rtmp3)*t8) if ((x*x+y*y).lt.potr*potr) then ir = ir+32 if (ir.gt.255) then ir = 255 endif ig = ig+32 if (ig.gt.255) then ig = 255 endif ib = ib+40 if (ib.gt.255) then ib = 255 endif ira2 = ira2+32.0d0*1.732d0 if (ira2.gt.255) then ira2 = 255 endif endif write(nfxv,"(3i4)") ir,ig,ib write(nf,"(3i4)") ira2 23022 continue 23020 continue close(unit = nf , status = "keep") close(unit = nfxv , status = "keep") call system('gzip fort.????') write(6,*) s endif ************************************************************ m0 = -m0 c*********< x. >******************************************* m1 = -m0 do 25024 j = 3,ly+2 m1 = -m1 if (m1.eq.1) then cpsie3 = cal*cpsi(3,j)+cbt*cpsi(2,j) i =2 icountmax = lxh else cpsie3 = cal*cpsi(2,j)+cbt*cpsi(1,j) i = 1 icountmax = lxh+1 endif do 25026 icount = 1,icountmax i = i+2 im = i-1 ip = i+1 cpsie1 = cpsie3 ctmp01 = hlf*(cpsi(i,j)+cpsi(ip,j)) ctmp02 = cgm*(cpsi(i,j)-cpsi(ip,j)) cpsie2 = ctmp01+ctmp02 cpsie3 = ctmp01-ctmp02 ctmp01 = hlf*(cpsie1+cpsie2) ctmp02 = cgm*(cpsie1-cpsie2) cpsi(im,j) = ctmp01+ctmp02 cpsi(i,j) = ctmp01-ctmp02 25026 continue 25024 continue c do 25110 j = 3,ly+2 cpsi(1,j) = cpsi(lx+1,j) cpsi(2,j) = cpsi(lx+2,j) cpsi(lx+3,j) = cpsi(3,j) cpsi(lx+4,j) = cpsi(4,j) 25110 continue do 27610 i = 3,lx+2 cpsi(i,1) = cpsi(i,ly+1) cpsi(i,2) = cpsi(i,ly+2) cpsi(i,ly+3) = cpsi(i,3) cpsi(i,ly+4) = cpsi(i,4) 27610 continue c***********< y. >******************************************* m1 = -m0 do 27024 i = 3,lx+2 m1 = -m1 if (m1.eq.1) then cpsie3 = cal*cpsi(i,3)+cbt*cpsi(i,2) j = 2 jcountmax = lyh else cpsie3 = cal*cpsi(i,2)+cbt*cpsi(i,1) j = 1 jcountmax = lyh+1 endif do 27026 jcount = 1,jcountmax j = j+2 jm =j-1 jp =j+1 cpsie1 = cpsie3 ctmp01 = hlf*(cpsi(i,j)+cpsi(i,jp)) ctmp02 = cgm*(cpsi(i,j)-cpsi(i,jp)) cpsie2 = ctmp01+ctmp02 cpsie3 = ctmp01-ctmp02 ctmp01 = hlf*(cpsie1+cpsie2) ctmp02 = cgm*(cpsie1-cpsie2) cpsi(i,jm) = ctmp01+ctmp02 cpsi(i,j) = ctmp01-ctmp02 27026 continue 27024 continue c do 27110 i = 3,lx+2 cpsi(i,1) = cpsi(i,ly+1) cpsi(i,2) = cpsi(i,ly+2) cpsi(i,ly+3) = cpsi(i,3) cpsi(i,ly+4) = cpsi(i,4) 27110 continue do 25510 j = 3,ly+2 cpsi(1,j) = cpsi(lx+1,j) cpsi(2,j) = cpsi(lx+2,j) cpsi(lx+3,j) = cpsi(3,j) cpsi(lx+4,j) = cpsi(4,j) 25510 continue c---------------< V. >---------------------------------- do 36423 j = 3,ly+2 do 36433 i = 3,lx+2 cpsi(i,j) = cev0(i,j)*cpsi(i,j) 36433 continue 36423 continue do 37880 i = 3,lx+2 cpsi(i,1) = cpsi(i,ly+1) cpsi(i,2) = cpsi(i,ly+2) cpsi(i,ly+3) = cpsi(i,3) cpsi(i,ly+4) = cpsi(i,4) 37880 continue do 35650 j = 3,ly+2 cpsi(1,j) = cpsi(lx+1,j) cpsi(2,j) = cpsi(lx+2,j) cpsi(lx+3,j) = cpsi(3,j) cpsi(lx+4,j) = cpsi(4,j) 35650 continue 53016 continue stop end