parameter (N=64, NN=2*N+1) implicit real*8(a-h,o-z) real x(0:NN), ubar(0:NN,0:3),uapp(0:NN) real uapp1(0:NN), xL(0:NN,0:3), ubars(0:NN,0:3), si(0:NN) complex abar(0:NN), al(0:NN) complex w, A, y, ss, ss1, ss2, F, sum, sum1 pi=4.0*atan(1.0) w=cmplx(0, 1) alpha=19.0*log(10.0) dx=2.0*pi/NN dx1=0.5*dx N2=0.0*N dt=0.01 ttime=0.0 Tmax=0.4 c compute the initial fourier coefficient do j=0, NN x(j)=j*dx if ( (x(j)-dx1) .lt. 0.9 .and. (x(j)+dx1) .gt. 0.9) jp=j enddo do j=0, NN s1=cos(0.5*(x(j)+dx1)) s2=cos(0.5*(x(j)-dx1)) if (j .lt. jp) then ubar(j, 0)=-2.0*(s1-s2)/dx elseif (j .gt. jp) then ubar(j, 0)=2.0*(s1-s2)/dx else p1=-2.0*(cos(0.45)-s2)/dx p2=2.0*(s1-cos(0.45))/dx ubar(j, 0)=p1+p2 endif s1=-cos(x(j)+dx1) s2=cos(x(j)-dx1) c ubar(j, 0)=0.3+0.7*(s1+s2)/dx enddo 99 do NR=0, 2 do k=0, NN-1 k1=k-0.5*(NN-1) ss=0.0 do j=0, NN-1 p1=k1*x(j) s1=cos(p1) s2=-sin(p1) ss=ss+ubar(j,NR)*cmplx(s1,s2) enddo abar(k)=ss/NN enddo do k=0, NN-1 k1=k-0.5*(NN-1) if (k1 .eq. 0) then si(k)=1.0 else p1=0.5*k1*dx p2=sin(p1) si(k)=p2/p1 endif al(k)=abar(k)/si(k) c ss1=-w*(k+0.5)*0.9 c ss1=exp(ss1) c ss2=-w*(k-0.5)*0.9 c ss2=exp(ss2) c ss1=(-ss1/(k1+0.5)+ss2/(k1-0.5)) c print *, al(k), ss1, NR, 'ppppppppppppp' enddo c compute the shock location N1=sqrt(1.0*N) N1=N1*sqrt(1.0*N1) N1=(1.0*N)**0.75 c N1=28 ss=(N1+1)*al(N+N1+1)/((N1+2)*al(N+N1+2)) y=-w*log(ss) ss=w*(N1+2)*y A=w*(N1+2)*al(N+N1+2)*exp(ss) if (NR .eq. 2) print *, 'haaaaaaaaaaa', A, y if (y1 .gt. 2.*pi .or. y1 .lt. 0.0) then y1=0.0 y=0.0 endif y1=real(y) y=real(y) A=real(A) do j=0, NN-1 xj=(j+0.5)*dx if (xj .le. y1 ) then F=-A*xj else F=A*(2*pi-xj) endif sum=0.0 do k=1, N if ( k .le. N2) then tau=1.0 else tau=exp(-alpha*(1.0*(k-N2)/N)**6) endif ps=k*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k+N)+w*A*ss1/k)*ss2 sum=sum+tau*ss k1=-k ps=k1*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(al(k1+N)+w*A*ss1/k1)*ss2 sum=sum+tau*ss enddo uapp1(j)=real(al(N)-A*(pi-y)+sum+F) c uapp1(j)=al(N)+sum+F enddo do j=0, NN if ( (x(j)-dx1) .lt. y1 .and. (x(j)+dx1) .gt. y1) ip=j enddo do j=0, NN sum1=0.0 xj=x(j) if (j .lt. ip) then F=-A*xj elseif ( j .gt. ip) then F=A*(2*pi-xj) elseif (j .eq. ip) then s1=y1**2-(xj-dx1)**2 s2=2.0*pi*(xj+dx1-y1) s3=(xj+dx1)**2-y1**2 F=A*(-0.5*s1+s2-0.5*s3)/dx endif do k=1, N if ( k .le. N2) then tau=1.0 else tau=exp(-alpha*(1.0*(k-N2)/N)**6) endif ps=k*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(abar(k+N)+w*si(k+N)*A*ss1/k)*ss2 sum1=sum1+tau*ss k1=-k ps=k1*y1 s1=cos(ps) s2=-sin(ps) ss1=cmplx(s1, s2) ps=k1*xj s1=cos(ps) s2=sin(ps) ss2=cmplx(s1, s2) ss=(abar(k1+N)+w*A*si(k1+N)*ss1/k1)*ss2 sum1=sum1+tau*ss enddo uapp(j)=real(abar(N)-A*(pi-y1)+sum1+F) enddo c Runge-Kutta do j=0, NN-1 ss1=uapp1(j) if (j .eq. 0) then ss2=uapp1(NN-1) else ss2=uapp1(j-1) endif xL(j, NR)=-0.5*(ss1*ss1-ss2*ss2)/dx c ubars(j, NR)=uapp(j) ubars(j, NR)=ubar(j, NR) if (NR .eq. 0) then ubar(j,1)=ubars(j,0)+dt*xL(j,0) endif if (NR .eq. 1) then sa=0.75*ubars(j,0)+0.25*ubars(j,1)+0.25*dt*xL(j,1) ubar(j,2)=sa endif if (NR .eq. 2) then sa=ubars(j,0)/3.0+ubars(j,2)*2/3.0+dt*xL(j,2)*2/3.0 ubar(j,3)=sa endif enddo enddo ttime=ttime+dt print *, ttime, ' time' if (ttime .lt. Tmax) then do j=0, NN ubar(j,0)=ubar(j, 3) enddo goto 99 endif write (17, *) 'x=[' do j=0, NN-1 write (17, *) x(j) enddo write (17, *) '];' write (18, *) 'y=[' do j=0, NN-1 write (18, *) real(ubar(j, 3)) enddo write (18, *) '];' print *, A, y, ' A, y' stop end