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 ENTERBOUNDARYVALUES(PHI,XSTART,DX,NX,YSTART,DY,NY)

!PRINTING OUT TO CONSOLE, TEMPORARY
COUNTER=1
DO WHILE(COUNTER<=NY+1)
WRITE(*,*) PHI((NY+2)-COUNTER,1:NX+1)
COUNTER=COUNTER+1
END DO
!ENDPRINT

CALL CREATEEQUATIONS(PHI, NX, NY, DX, DY, EQS, XSTART, YSTART)
WRITE(*,*) 'Har kommet hit også gitt'
END PROGRAM


!---------------------------------!
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
WRITE(*,*) '1'
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
WRITE(*,*) '2'
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
WRITE(*,*) '3'
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
WRITE(*,*) '4'

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 TRV
REAL, DIMENSION(:,:), ALLOCATABLE :: TMP
REAL, DIMENSION(:,:), ALLOCATABLE :: FACTORS

INTEGER TMPI
ALLOCATE (TMP(1,((NX-1)*(NY-1)+1)))
ALLOCATE (FACTORS((NX-1)*(NY-1)))
DO TRV=1, ((NY-1)*(NX-1))
FACTORS=0.0
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