PROGRAM MYPROG
INTEGER NX, NY ! Variables to hold number of intervals in each direction
REAL XSTART, XEND ! Variables defining the interval in x-direction
REAL YSTART, YEND ! Variables defining the inverval in y-direction
REAL DX, DY
!REAL PHI(21,21)
!REAL EQS((20*20),(20*20)+1)
INTEGER COUNTER ! TEMPORARY, USED FOR PRINTING VALUES
REAL, DIMENSION(: , : ), ALLOCATABLE :: PHI
REAL, DIMENSION(: , : ), ALLOCATABLE :: EQS
CALL DEFINEMATRIX(XSTART, XEND, YSTART, YEND, NX, NY) ! Lets the user specify the interval for which the problem will be solved, and the number of intervals for discretisation
ALLOCATE (PHI(NY+1,NX+1))
ALLOCATE (EQS((NY-1)*(NX-1),((NY-1)*(NX-1))+1))
DX=(XEND-XSTART)/NX
DY=(YEND-YSTART)/NY
CALL PRINT(PHI,NX,NY)
WRITE(*,*) 'Før verdier er ført inn'
CALL ENTERBOUNDARYVALUES(PHI,XSTART,DX,NX,YSTART,DY,NY)
CALL PRINT(PHI,NX,NY)
WRITE(*,*) 'Med grenseverdier'
CALL CREATEEQUATIONS(PHI, NX, NY, DX, DY, EQS, XSTART, YSTART)
WRITE(*,*) 'Har kommet hit også gitt'
CALL GAUSS(EQS, NX, NY)
CALL ENTERINNERVALUES(PHI, NX, NY, EQS)
CALL PRINT(PHI,NX,NY)
WRITE(*,*) 'Med indre verdier'
CALL CREATEOUTPUT(PHI, NX, NY)
WRITE(*,*) 'Verdier skrevet til output.dat'
END PROGRAM
SUBROUTINE PRINT(MATRIX, NX, NY)
REAL MATRIX(NY+1,NX+1)
INTEGER NX, NY
INTEGER COUNTER
!PRINTING OUT TO CONSOLE, TEMPORARY
COUNTER=1
DO WHILE(COUNTER<=NY+1)
WRITE(*,*) MATRIX((NY+2)-COUNTER,1:NX+1)
COUNTER=COUNTER+1
END DO
!ENDPRINT
END SUBROUTINE
!---------------------------------!
SUBROUTINE DEFINEMATRIX(XS,XE,YS,YE,NX,NY)
REAL XS, XE, YS, YE
INTEGER NX, NY
WRITE(*,*) 'Choose the interval in x-direction for which you wish to solve the problem. Input format <x_min> <x_max>'
CALL READINTERVAL(XS, XE) ! Prompt user to input x-interval
WRITE(*,*) 'Choose the interval in y-direction for which you wish to solve the problem. Input format <y_min> <y_max>'
CALL READINTERVAL(YS, YE) ! Prompt user to input y-interval
WRITE(*,*) 'Specify number of intervals you want in the x-direction'
READ(*,*) NX
WRITE(*,*) 'Specify number of intervals you want in the y-direction'
READ(*,*) NY
END SUBROUTINE
!----------------------------------!
SUBROUTINE READINTERVAL(START, END)
REAL START, END
read(*,*) START, END ! Assignes the values to the corresponding variables
END SUBROUTINE
!----------------------------------!
!----------------------------------!
REAL FUNCTION PHIBOUNDARY(X,Y)
REAL X, Y
PHIBOUNDARY=(X**2+Y**2)/4.0
RETURN
END FUNCTION
!----------------------------------!
SUBROUTINE ENTERBOUNDARYVALUES(PHI, XSTART, DX, NX, YSTART, DY, NY)
REAL PHI(NY+1,NX+1)
REAL XSTART, DX, YSTART, DY
INTEGER NX, NY
INTEGER COL
INTEGER ROW
DO WHILE(COL<=NX+1)
PHI(1,COL)=PHIBOUNDARY(XSTART+(COL-1)*DX,YSTART)
PHI(NY+1,COL)=PHIBOUNDARY(XSTART+(COL-1)*DX,YSTART+NY*DY)
COL=COL+1
END DO
ROW=2
DO WHILE(ROW<=NY)!DENNE GÅR FRA 2 TIL OG MED NY, DA 1 OG NY+1 ER FYLT UT AV LØKKEN OVER
PHI(ROW,1)=PHIBOUNDARY(XSTART,YSTART+(ROW-1)*DY)
PHI(ROW,NX+1)=PHIBOUNDARY(XSTART+NX*DX,YSTART+(ROW-1)*DY)
ROW=ROW+1
END DO
END SUBROUTINE
REAL FUNCTION GET1DPOS(COL,ROW,NX)
INTEGER COL, ROW, NX
GET1DPOS=INT(INT(COL-1)+INT(INT(ROW-2)*INT(NX-1)))
RETURN
END FUNCTION
REAL FUNCTION G(X,Y)
REAL X,Y
G=1.0
RETURN
END FUNCTION
SUBROUTINE CREATEEQUATIONS(PHI, NX, NY, DX, DY, EQS, XSTART, YSTART)
REAL PHI(NY+1,NX+1)
INTEGER NX, NY
REAL DX, DY
REAL EQS((NX-1)*(NY-1),((NX-1)*(NY-1))+1)
REAL XSTART, YSTART
REAL, DIMENSION(: , : ), ALLOCATABLE :: EQ
INTEGER POSY, POSX
INTEGER ROW, COL
REAL COEFFDY, COEFFDX, COEFFG
REAL STATICCOEFF
INTEGER IND
ALLOCATE (EQ(1,((NX-1)*(NY-1))+1))
COEFFDY=-1.0*(DX**2)/(2.0*((DX**2)+(DY**2)))
COEFFDX=-1.0*(DY**2)/(2.0*((DX**2)+(DY**2)))
COEFFG=(DX**2)*(DY**2)/(2.0*((DX**2)+(DY**2)))
EQS=0.0
POSY=2
DO WHILE(POSY<=NY)
POSX=2
DO WHILE(POSX<=NX)
EQ=0.0
STATICCOEFF=0.0
WRITE(*,*) 'Here we go once more (',POSY,',',POSX,')'
IF (POSY.EQ.2) THEN
STATICCOEFF=STATICCOEFF-1.0*COEFFDY*PHI(POSY-1,POSX)
ELSE
IND=INT(GET1DPOS(POSX,POSY-1,NX))
EQ(1,IND)=COEFFDY
END IF
IF (POSX.EQ.2) THEN
STATICCOEFF=STATICCOEFF-1.0*COEFFDX*PHI(POSY,POSX-1)
ELSE
IND=INT(GET1DPOS(POSX-1,POSY,NX))
EQ(1,IND)=COEFFDX
END IF
IF (POSY.EQ.NY) THEN
STATICCOEFF=STATICCOEFF-1.0*COEFFDY*PHI(POSY+1,POSX)
ELSE
IND=INT(GET1DPOS(POSX,POSY+1,NX))
EQ(1,IND)=COEFFDY
END IF
IF (POSX.EQ.NX) THEN
STATICCOEFF=STATICCOEFF-1.0*COEFFDX*PHI(POSY,POSX+1)
ELSE
IND=INT(GET1DPOS(POSX+1,POSY,NX))
EQ(1,IND)=COEFFDX
END IF
STATICCOEFF=STATICCOEFF+COEFFG*G((XSTART+(POSX-1)*DX),(YSTART+(POSY-1)*DY))
EQ(1,((NX-1)*(NY-1))+1)=STATICCOEFF
IND=INT(GET1DPOS(POSX,POSY,NX))
EQ(1,IND)=1
EQS(IND,:)=EQ(1,:)
WRITE(*,*) EQ(1,:)
POSX=POSX+1
END DO
POSY=POSY+1
END DO
!TEMP
IND=1
WRITE(*,*) 'Dette er en test av lignignsmatrisa'
DO WHILE(IND<=(NX-1)*(NY-1))
WRITE(*,*) EQS(IND,:)
IND=IND+1
END DO
WRITE(*,*) 'Reached the end, whaddya know'
!ENDTEMP
END SUBROUTINE
SUBROUTINE GAUSS(EQS, NX, NY)
REAL EQS((NX-1)*(NY-1),((NX-1)*(NY-1)+1))
INTEGER NX, NY
INTEGER ROW, COL
INTEGER RC
REAL, DIMENSION(:,:), ALLOCATABLE :: TMP
REAL FACTOR
INTEGER IND
!REAL, DIMENSION(:,:), ALLOCATABLE :: FACTORS
INTEGER TMPI
ALLOCATE (TMP(1,((NX-1)*(NY-1)+1)))
!ALLOCATE (FACTORS((NX-1)*(NY-1)))
WRITE(*,*) 'Gausselimination commences!'
DO RC=1, ((NY-1)*(NX-1))
DO ROW=1, ((NY-1)*(NX-1))
IF (EQS(RC,RC).EQ.0) THEN
WRITE(*,*) '0 på RC, bytt om rader'
END IF
IF (.NOT.(RC.EQ.ROW)) THEN
FACTOR=EQS(ROW,RC)/EQS(RC,RC)
!DO COL=RC, ((NY-1)*(NX-1)+1)
!EQS(ROW,COL)=EQS(ROW,COL)-FACTOR*EQS(RC,COL)
!END DO
EQS(ROW,:)=EQS(ROW,:)-FACTOR*EQS(RC,:)
END IF
END DO
END DO
!TEMP
IND=1
WRITE(*,*) 'Dette er en test av gausseliminert'
DO WHILE(IND<=(NX-1)*(NY-1))
WRITE(*,*) EQS(IND,:)
IND=IND+1
END DO
WRITE(*,*) 'Reached the end, whaddya know'
!ENDTEMP
!DO ROW=1, ((NY-1)*(NX-1))
!IF (ROW==TRV) THEN
!FACTORS(ROW)=0.0
!ELSE
!FACTORS(ROW)=EQS(ROW,TRV)/EQS(TRV,TRV)
!DO COL=1, ((NY-1)*(NX-1))
!DO ROW=1, ((NY-1)*(NX-1))
END SUBROUTINE
SUBROUTINE ENTERINNERVALUES(PHI,NX,NY,EQS)
REAL PHI(NY+1,NX+1)
INTEGER NX, NY
REAL EQS((NY-1)*(NX-1),((NY-1)*(NX-1))+1)
INTEGER RC
INTEGER MODU
INTEGER ROWP, COLP
INTEGER NUMROW
DO RC=1, ((NY-1)*(NX-1))
NUMROW=INT((RC-1)/(NX-1))
ROWP=2+NUMROW
COLP=1+RC-NUMROW*(NX-1)
PHI(ROWP,COLP)=(EQS(RC,((NY-1)*(NX-1)+1))/EQS(RC,RC))
WRITE(*,*) 'Row',ROWP,'Col',COLP,'Val',(EQS(RC,((NY-1)*(NX-1)+1))/EQS(RC,RC))
!MODU=MOD((RC-1),(NX-1))
!COLP=2+MODU
!ROWP=2-INT((RC-1)/(NX-1))
!PHI(ROWP,COLP)=EQS(RC,((NY-1)*(NX-1)+1))/EQS(RC,RC)
!MODU=MOD((RC-1),(NX-1))
!ROWP=2+MODU
!COLP=1+RC-(MODU*(NX-1))
END DO
END SUBROUTINE
SUBROUTINE CREATEOUTPUT(PHI, NX, NY)
REAL PHI(NY+1,NX+1)
INTEGER NX
INTEGER NY
INTEGER COUNTX
INTEGER COUNTY
OPEN(UNIT=11, FILE='output.dat', ACTION='write', STATUS='unknown')
WRITE(11,*), NX, NY
DO COUNTY=1, NY+1
DO COUNTX=1, NX+1
WRITE(11,*) PHI(COUNTY,COUNTX)
END DO
END DO
CLOSE(11)
!" nx ny
! phi(0*dx,0)
! phi(1*dx,0)
! ......
! phi(nx*dx,0)
! phi(0*dx,dx)
! phi(1*dx,dx)
! .........
!"
END SUBROUTINE