awp.f90

PROGRAM awp

IMPLICIT NONE

INTEGER:: N, i, j, idat
DOUBLE PRECISION:: t0, tend, E_1, E_2, dt_1, dt_2, o,dy(4),y(4),t,yex(2),v0,r0,T_0,pi
CHARACTER::datei
pi=acos(-1.)

t0=0.
tend=0.6
v0=10.
r0=1.
T_0=2.*pi*r0*1./v0


OPEN(UNIT=100,FILE='ausgabe_ec.csv',STATUS='UNKNOWN')
WRITE(100,'(A)') "Streifen; analytisch; numerisch; Schrittweite; Fehler; Ordnung"

idat=1
DO i=4,10
    N=2**i
    dt_2=DBLE(tend-t0)/DBLE(N)
    y=(/ r0,0.,0.,v0 /)
    t=0.
    E_2=0.

    WRITE(datei,'(I1)')idat
    idat=idat+1
    write(*,*)idat
    OPEN(UNIT=101,FILE='ec_erg'//TRIM(datei)//'.csv',STATUS='UNKNOWN')
    !OPEN(UNIT=101,FILE='rk_erg'//TRIM(datei)//'.csv',STATUS='UNKNOWN')
    WRITE(101,*)"Zeit;x;xex;y;yex;vx;vy"

    !Anfangswerte
    yex=yexakt(t,r0,T_0)
    WRITE(101,'(F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8)') &
      & t,";",y(1),";",yex(1),";",y(2),";",yex(2),";",y(3),";",y(4)

    DO j=0,N-1
        dy = eulercauchy(y,dt_2,t)
        !dy = rungekutta(y,dt,t)
        y = y + dy
        t = t + dt_2
        yex=yexakt(t,r0,T_0)
        WRITE(101,'(F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8)') &
          & t,";",y(1),";",yex(1),";",y(2),";",yex(2),";",y(3),";",y(4)

        E_2 = E_2 + (y(1)-yex(1))**2
    END DO
    CLOSE(101)

    E_2=sqrt(E_2/N)

    IF (i>4) THEN
        o=DLOG(E_1/E_2)/DLOG(dt_1/dt_2)
    ELSE
        o=0.
    END IF

    WRITE(100,'(I8, A, F12.8, A, F12.8, A, F12.8, A, E14.7, A, F10.3)')N,&
    &";", yex(1), ";", y(1), ";", dt_2, ";", E_2, ";", o

    E_1=E_2
    dt_1=dt_2
END DO

CLOSE(100)

write(*,*) y(1)
write(*,*) y(2)
write(*,*) y(3)
write(*,*) y(4)

CONTAINS





function f(y,t)
    real::f(4)
    real::B
    real, intent(in)::y(4),t
    B=-10.
    write(*,*) "vx,vy: ",y(3)," ",y(4)
    f=(/ y(3),y(4),y(4)*B,-y(3)*B /)
end function f

function yexakt(t,r,T_0)
    real::yexakt(2),pi
    real, intent(in)::t,r,T_0
    pi=acos(-1.)
    yexakt=(/ r*cos(2.*pi/T_0*t), r*sin(2.*pi/T_0*t) /)
end function yexakt


function eulercauchy(y,dt,t)
    real::eulercauchy(4)
    real, intent(in) :: y(4), dt, t
    eulercauchy=dt*f(y,t)
end function eulercauchy

function rungekutta(y,dt,t)
    real::rungekutta(4)
    real, intent(in) :: y(4), dt, t
    real:: k1(4), k2(4), k3(4), k4(4)

    k1=f(y,t)
    k2=f(y+dt/2*k1,t+dt/2.)
    k3=f(y+dt/2*k2,t+dt/2.)
    k4=f(y+dt*k3,t+dt)
    rungekutta=dt/6.*(k1+2.*(k2+k3)+k4)

end function rungekutta


END PROGRAM