module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defpi


module cdraw_data
    implicit none
    integer::ind,ibc,ifc,isc
    real(4)::baselength,grdspc
    real(4)::xmin,xmax,ymin,ymax
    real(4)::bav,bfv
    real(4)::conint,annint
    character::strx*100,stry*100
    integer::ncut,numcut(10)

    real(4)::fmax,fact0=1.0E-10
end module cdraw_data


program f90_gmt_ini4_vec
    !Drawing by GMT for 2D Stress Analysis
    use defpi
    use cdraw_data
    implicit none
    integer::i,j,k,l,kk,ne,n,ia,ja,k1,k2,n1,n2,n3
    character::strcom*50                !Comment
    integer::NSTRES                     !Stress state (0: plane strain, 1: plane stress)
    integer::NODT                       !Number of nodes
    integer::NELT                       !Number of elements
    integer::MATEL                      !Number of material sets
    integer::KOX                        !Number of restricted nodes in x-direction
    integer::KOY                        !Number of restricted nodes in y-direction
    integer::KXY                        !Number of restricted nodes in x and y-directions
    integer::NF                         !Number of loaded nodes
    integer::IPR
    integer,allocatable::kakom(:,:)     !Element connectivity
    integer,allocatable::matno(:)       !Material set number
    real(4),allocatable::x(:)           !x-coordinate of node
    real(4),allocatable::y(:)           !y-coordinate of node
    integer,allocatable::nokx(:)        !Restricted node number in x-direction
    integer,allocatable::noky(:)        !Restricted node number in y-direction
    integer,allocatable::noxy(:)        !Restricted node number in x and y-directions
    real(4),allocatable::f(:)           !External load vector
    real(4),allocatable::dis(:)         !Displacement vector
    real(4),allocatable::strsave(:,:,:) !Stresses of element
    integer::nt                         !Total degree of freedom (NODT x nhen)
    integer::nod                        !Number of nodes per element
    integer::nhen=2                     !Number of degree of freedom per node
    integer::kkmax                      !Number of integration point
    real(4),allocatable::xg(:,:)        !x-coordinate of centroid of element
    real(4),allocatable::yg(:,:)        !y-coordinate of centroid of element
    real(4),allocatable::sf(:,:)        !Safety Factor
    real(4)::rw
    integer::iw

    character::linebuf*1000
    character::fnameR*50

    call getarg(1,fnameR) !Input file name

    !***************************
    !data input
    !***************************
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,'()')
    read(11,*) nod,NODT,NELT,MATEL,KOX,KOY,NF,NSTRES,IPR
    nt=NODT*nhen
    kkmax=4
    if(IPR==1)kkmax=1
    allocate(kakom(1:NELT,1:nod))
    allocate(matno(1:NELT))
    allocate(x(1:NODT))
    allocate(y(1:NODT))
    allocate(nokx(1:KOX))
    allocate(noky(1:KOY))
    allocate(noxy(1:KOX))
    allocate(f(1:nt))
    allocate(dis(1:nt))
    allocate(strsave(1:NELT,1:kkmax,1:3))
    allocate(xg(1:NELT,1:kkmax))
    allocate(yg(1:NELT,1:kkmax))
    allocate(sf(1:NELT,1:kkmax))
    read(11,'()')
    read(11,'()')
    n1=0
    n2=0
    n3=0
    do i=1,NODT
        read(11,*) iw,rw,rw,f(2*i-1),f(2*i),k1,k2,rw,rw
        if(0<KOX)then
            if(k1==1.and.k2==0)then
                n1=n1+1
                nokx(n1)=i
            end if
        end if
        if(0<KOY)then
            if(k1==0.and.k2==1)then
                n2=n2+1
                noky(n2)=i
            end if
        end if
        if(0<KOX.and.0<KOY)then
            if(k1==1.and.k2==1)then
                n3=n3+1
                noxy(n3)=i
            end if
        end if
    end do
    KOX=n1
    KOY=n2
    KXY=n3
    read(11,'()')
    read(11,'()')
    do ne=1,NELT
        read(11,*) iw,(kakom(ne,j),j=1,nod)
    end do
    read(11,'()')
    read(11,'()')
    do i=1,NODT
        read(11,*) iw,x(i),y(i),dis(2*i-1),dis(2*i),rw,rw,rw,rw
    end do
    read(11,'()')
    read(11,'()')
    do ne=1,NELT
        do kk=1,kkmax
            read(11,*) iw,iw,xg(ne,kk),yg(ne,kk),(strsave(ne,kk,j),j=1,3),rw,rw,rw,sf(ne,kk),matno(ne)
        end do
    end do
    close(11)

    open(11,file='_inp_base.txt',status='old')
    read(11,*) ind,ibc,ifc,isc
    read(11,*) baselength,grdspc
    read(11,*) bav,bfv
    read(11,*) conint,annint
    read(11,*) strx
    read(11,*) stry
    read(11,*) ncut,(numcut(i),i=1,ncut)
    close(11)

    xmin=(x(1))
    xmax=(x(1))
    ymin=(y(1))
    ymax=(y(1))
    do i=1,NODT
        if((x(i))<xmin)xmin=(x(i))
        if(xmax<(x(i)))xmax=(x(i))
        if((y(i))<ymin)ymin=(y(i))
        if(ymax<(y(i)))ymax=(y(i))
    end do

    if(0<ind)call DAT_IND(kakom,x,y,NELT,NODT,nod)
    call DAT_MESH(kakom,matno,x,y,nokx,noky,noxy,f,nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt)
    call DAT_VEC(xg,yg,strsave,NELT,kkmax)
    call MAKE_BAT(MATEL,KOX,KOY,KXY)
    call system('bash _gmt_vec.txt')

stop

contains

    subroutine DAT_IND(kakom,x,y,NELT,NODT,nod)
        integer,intent(in)::NELT,NODT,nod
        integer,intent(in)::kakom(1:NELT,1:nod)
        real(4),intent(in)::x(1:NODT),y(1:NODT)

        integer::ne,i,j,k,l
        real(4)::xg,yg,dts

        dts=0.001*baselength

        open (12, file='_ne_circ.txt', status='replace')
        open (13, file='_text_ne.txt', status='replace')
        do ne=1,NELT
            i=kakom(ne,1)
            j=kakom(ne,2)
            k=kakom(ne,3)
            l=kakom(ne,4)
            xg=0.25*(x(i)+x(j)+x(k)+x(l))
            yg=0.25*(y(i)+y(j)+y(k)+y(l))
            write(12,*) xg,' ',yg
            write(13,*) xg,' ',yg,' 9 0 MC ',ne
        end do
        close(12)
        close(13)

        open (12, file='_text_nd.txt', status='replace')
        do i=1,NODT
            write(12,*) x(i)-dts,' ',y(i)-dts,' 9 0 TR ',i
        end do
        close(12)
    end subroutine DAT_IND


    subroutine DAT_MESH(kakom,matno,x,y,nokx,noky,noxy,f,nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt)
        integer::nod,NODT,NELT,MATEL,KOX,KOY,KXY,nt
        integer::kakom(1:NELT,1:nod),matno(1:NELT)
        integer::nokx(1:KOX),noky(1:KOY),noxy(1:KOX)
        real(4)::x(1:NODT),y(1:NODT),f(1:nt)

        integer::i,j,k,l,ma,ne
        real(4)::fscale,dscale,dismax
        character::filename*50

        if(ymax-ymin<xmax-xmin)then
            fscale=0.05*(xmax-xmin)
        else
            fscale=0.05*(ymax-ymin)
        end if

        fmax=0.0
        do i=1,NODT
            if(fmax<sqrt(f(2*i-1)**2.0+f(2*i)**2.0))fmax=sqrt(f(2*i-1)**2.0+f(2*i)**2.0)
        end do

        do ma=1,MATEL
            write (filename, '("_mesh", i2.2, ".txt")') ma
            open (12, file=filename, status='replace')
            do ne=1,NELT
                if(matno(ne)==ma)then
                    i=kakom(ne,1)
                    j=kakom(ne,2)
                    k=kakom(ne,3)
                    l=kakom(ne,4)
                    write(12,'(a)') '>'
                    write(12,*) x(i),' ',y(i)
                    write(12,*) x(j),' ',y(j)
                    write(12,*) x(k),' ',y(k)
                    write(12,*) x(l),' ',y(l)
                end if
            end do
            close (12)
        end do

        if(0<KOX)then
            open (12, file='_kox.txt', status='replace')
            do i=1,KOX
                write(12,*) x(nokx(i)),' ',y(nokx(i))
            end do
            close(12)
        end if
        if(0<KOY)then
            open (12, file='_koy.txt', status='replace')
            do i=1,KOY
                write(12,*) x(noky(i)),' ',y(noky(i))
            end do
            close(12)
        end if
        if(0<KXY)then
            open (12, file='_kxy.txt', status='replace')
            do i=1,KXY
                write(12,*) x(noxy(i)),' ',y(noxy(i))
            end do
            close(12)
        end if

        if(fact0<fmax)then
            open (12, file='_force.txt', status='replace')
            do i=1,NODT
                if(fact0<sqrt(f(2*i-1)**2.0+f(2*i)**2.0))then
                    write(12,*) x(i),' ',y(i),' ',x(i)+f(2*i-1)/fmax*fscale,' ',y(i)+f(2*i)/fmax*fscale
                end if
            end do
            close(12)
        end if

        dismax=0.0
        do i=1,NODT
            if(dismax<abs(dis(i)))dismax=abs(dis(i))
        end do
        open(12,file='_dis.txt',status='replace')
        dscale=0.01*(xmax-xmin)/dismax
        do ne=1,NELT
            i=kakom(ne,1)
            j=kakom(ne,2)
            k=kakom(ne,3)
            l=kakom(ne,4)
            write(12,'(a)') '>'
            write(12,*) x(i)+dscale*dis(2*i-1),' ',y(i)+dscale*dis(2*i)
            write(12,*) x(j)+dscale*dis(2*j-1),' ',y(j)+dscale*dis(2*j)
            write(12,*) x(k)+dscale*dis(2*k-1),' ',y(k)+dscale*dis(2*k)
            write(12,*) x(l)+dscale*dis(2*l-1),' ',y(l)+dscale*dis(2*l)
        end do
        close(12)
    end subroutine DAT_MESH


    subroutine DAT_VEC(xg,yg,strsave,NELT,kkmax)
        integer::NELT,kkmax
        real(4)::xg(1:NELT,1:kkmax),yg(1:NELT,1:kkmax)
        real(4)::strsave(1:NELT,1:kkmax,1:3)
        integer::ne,kk
        real(4)::sigx,sigy,tauxy,ps1,ps2,ang
        real(4)::prst1(1:NELT),prst2(1:NELT),angle(1:NELT)
        real(4)::ps1max,ps2max,scl,dsmax
        real(4)::x1,y1,x2,y2
        real(4)::xxg(1:NELT),yyg(1:NELT)

        ps1max=0.0
        ps2max=0.0
        do ne=1,NELT
            sigx=0.0
            sigy=0.0
            tauxy=0.0
            xxg(ne)=0.0
            yyg(ne)=0.0
            do kk=1,kkmax
                sigx=sigx+strsave(ne,kk,1)
                sigy=sigy+strsave(ne,kk,2)
                tauxy=tauxy+strsave(ne,kk,3)
                xxg(ne)=xxg(ne)+xg(ne,kk)
                yyg(ne)=yyg(ne)+yg(ne,kk)
            end do
            sigx=sigx/real(kkmax)
            sigy=sigy/real(kkmax)
            tauxy=tauxy/real(kkmax)
            xxg(ne)=xxg(ne)/real(kkmax)
            yyg(ne)=yyg(ne)/real(kkmax)
            call PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
            prst1(ne)=ps1
            prst2(ne)=ps2
            angle(ne)=ang
            if(ps1max<abs(ps1))ps1max=abs(ps1)
            if(ps2max<abs(ps2))ps2max=abs(ps2)
        end do
        dsmax=max(xmax-xmin,ymax-ymin)
        open (12, file='_vec1b.txt', status='replace')
        open (13, file='_vec1r.txt', status='replace')
        scl=0.05*dsmax/abs(ps1max)
        do ne=1,ne
            ang=angle(ne)
            x1=xxg(ne)-0.5*prst1(ne)*cos(ang/180.0*pi)*scl
            x2=xxg(ne)+0.5*prst1(ne)*cos(ang/180.0*pi)*scl
            y1=yyg(ne)-0.5*prst1(ne)*sin(ang/180.0*pi)*scl
            y2=yyg(ne)+0.5*prst1(ne)*sin(ang/180.0*pi)*scl
            if(0.0<=prst1(ne))write(12,*) x1,' ',y1,' ',x2,' ',y2
            if(prst1(ne)<0.0) write(13,*) x1,' ',y1,' ',x2,' ',y2
        end do
        close(12)
        close(13)
        open (12, file='_vec2b.txt', status='replace')
        open (13, file='_vec2r.txt', status='replace')
        scl=0.05*dsmax/abs(ps2max)
        do ne=1,ne
            ang=angle(ne)
            x1=xxg(ne)-0.5*prst2(ne)*sin(ang/180.0*pi)*scl
            x2=xxg(ne)+0.5*prst2(ne)*sin(ang/180.0*pi)*scl
            y1=yyg(ne)+0.5*prst2(ne)*cos(ang/180.0*pi)*scl
            y2=yyg(ne)-0.5*prst2(ne)*cos(ang/180.0*pi)*scl
            if(0.0<=prst2(ne))write(12,*) x1,' ',y1,' ',x2,' ',y2
            if(prst2(ne)<0.0) write(13,*) x1,' ',y1,' ',x2,' ',y2
        end do
        close(12)
        close(13)

    end subroutine DAT_VEC


    subroutine MAKE_BAT(MATEL,KOX,KOY,KXY)
        integer::MATEL,KOX,KOY,KXY

        integer::ma
        character::srange*100,sscale*100,sba*100,sgrdspc*100,ssradi*100,filename*100
        character::ssx*100,ssy*100,scc*50,sca*50,sscld*100,scol*50,snum*50,spsize*50,scsize*50
        character::strw1*50,strw2*50,strw3*50,strw4*50
        real(4)::psize,ds
        real(4)::work
        character(len=3)::col_r(1:MATEL),col_g(1:MATEL),col_b(1:MATEL)
        character(len=11)::col_mesh(1:MATEL)

        psize=0.015*baselength
        ds=bav

        open (11, file='_col_mesh.txt', status='old')
        do i=1,MATEL
            read(11,*) col_r(i),col_g(i),col_b(i)
            col_mesh(i)=trim(adjustl(col_r(i)))//'/'//trim(adjustl(col_g(i)))//'/'//trim(adjustl(col_b(i)))
        end do
        close(11)

        write(strw1,'(f10.3)') xmin-ds;
        write(strw2,'(f10.3)') xmax;
        write(strw3,'(f10.3)') ymin-ds;
        write(strw4,'(f10.3)') ymax;
        srange='range='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/'//trim(adjustl(strw4))
        if(ymax-ymin<=xmax-xmin)then
            write(strw1,'(f10.3)') baselength;
            write(strw2,'(f10.3)') baselength/(xmax-xmin+ds)*(ymax-ymin+ds)
        else
            write(strw1,'(f10.3)') baselength/(ymax-ymin+ds)*(xmax-xmin+ds);
            write(strw2,'(f10.3)') baselength;
        end if
        sscale='scale='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') bav
        write(strw2,'(f10.3)') bfv
        sba='ba=a'//trim(adjustl(strw1))//'f'//trim(adjustl(strw2))
        write(strw1,'(f10.3)') grdspc
        sgrdspc='grdspc='//trim(adjustl(strw1))
        write(strw1,'(f11.4)') grdspc/sqrt(2.0)
        ssradi='sradi='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') conint
        scc='conint='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') annint
        sca='annint='//trim(adjustl(strw1))
        ssx='strx="'//trim(adjustl(strx))//'"'
        ssy='stry="'//trim(adjustl(stry))//'"'
        if(ymax-ymin<=xmax-xmin)then
            work=baselength/(xmax-xmin)*(ymax-ymin)
            write(strw1,'(f10.3)') baselength+0.5
            write(strw2,'(f10.3)') work-work*0.5/2.0-0.3
            write(strw3,'(f10.3)') work*0.5
            if(4.0<work*0.5)then
                write(strw2,'(f10.3)') work-4.0/2.0-0.3
                write(strw3,'(f10.3)') 4.0
            end if
        else
            work=baselength/(ymax-ymin+ds)*(xmax-xmin+ds)
            write(strw1,'(f10.3)') work+0.5
            write(strw2,'(f10.3)') baselength-work*0.5/2.0-0.3
            write(strw3,'(f10.3)') work*0.5
            if(4.0<work*0.5)then
                write(strw2,'(f10.3)') work-4.0/2.0-0.3
                write(strw3,'(f10.3)') 4.0
            end if
        end if
        sscld='scld='//trim(adjustl(strw1))//'/'//trim(adjustl(strw2))//'/'//trim(adjustl(strw3))//'/0.3'
        write(strw1,'(f10.3)') psize
        spsize='psize='//trim(adjustl(strw1))
        write(strw1,'(f10.3)') psize*4.0
        scsize='csize='//trim(adjustl(strw1))

        !writing of batch file
        open (12, file='_gmt_vec.txt', status='replace')
        write(12,'(a)') 'gmt set FONT_ANNOT_PRIMARY 10'
        write(12,'(a)') 'gmt set MAP_ANNOT_OFFSET_PRIMARY 0.3c'
        write(12,'(a)') 'gmt set FONT_LABEL 10p'
        write(12,'(a)') 'gmt set MAP_LABEL_OFFSET 0.3c'
        write(12,'(a)') 'gmt set MAP_TICK_LENGTH_PRIMARY 0.2c'
        write(12,'(a)') trim(adjustl(srange))
        write(12,'(a)') trim(adjustl(sscale))
        write(12,'(a)') trim(adjustl(sba))
        write(12,'(a)') trim(adjustl(sgrdspc))
        write(12,'(a)') trim(adjustl(ssradi))
        write(12,'(a)') trim(adjustl(scc))
        write(12,'(a)') trim(adjustl(sca))
        write(12,'(a)') trim(adjustl(ssx))
        write(12,'(a)') trim(adjustl(ssy))
        write(12,'(a)') trim(adjustl(sscld))
        write(12,'(a)') trim(adjustl(spsize))
        write(12,'(a)') trim(adjustl(scsize))
        write(12,'()')
        !drawing of ps1 vector
        write(12,'(a)') '# *** Drawing of vector ***'
        write(12,'(a)') 'fig=_fig_gmt_vec1.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$stry" -BWS -P -K > $fig'
        do ma=1,MATEL
            write (strw1,'(i2.2)') ma
            snum='num='//trim(adjustl(strw1))
            scol='col='//col_mesh(ma)
            write(12,'(a)') trim(adjustl(snum))
            write(12,'(a)') trim(adjustl(scol))
            write(12,'(a)') 'gmt psxy _mesh${num}.txt -R -J -W0.5,#aaaaaa -G$col -K -O >> $fig'
        end do
        write(12,'(a)') 'gmt psxy _vec1b.txt -R -J -Sv0.2+s -W0.5,blue -K -O >> $fig'
        write(12,'(a)') 'gmt psxy _vec1r.txt -R -J -Sv0.2+s -W0.5,red  -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        !drawing of ps2 vector
        write(12,'(a)') 'fig=_fig_gmt_vec2.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$stry" -BWS -P -K > $fig'
        do ma=1,MATEL
            write (strw1,'(i2.2)') ma
            snum='num='//trim(adjustl(strw1))
            scol='col='//col_mesh(ma)
            write(12,'(a)') trim(adjustl(snum))
            write(12,'(a)') trim(adjustl(scol))
            write(12,'(a)') 'gmt psxy _mesh${num}.txt -R -J -W0.5,#aaaaaa -G$col -K -O >> $fig'
        end do
        write(12,'(a)') 'gmt psxy _vec2b.txt -R -J -Sv0.2+s -W0.5,blue -K -O >> $fig'
        write(12,'(a)') 'gmt psxy _vec2r.txt -R -J -Sv0.2+s -W0.5,red  -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
        !drawing of displacement mode
        write(12,'(a)') '# *** Drawing of displacement mode ***'
        write(12,'(a)') 'range=-10/10/-10/10'
        write(12,'(a)') 'fig=_fig_gmt_disp.eps'
        write(12,'(a)') 'gmt psbasemap -R$range -JX$scale -Bx$ba+l"$strx" -By$ba+l"$stry" -BWS -P -K > $fig'
        do ma=1,MATEL
            write (strw1,'(i2.2)') ma
            snum='num='//trim(adjustl(strw1))
            scol='col='//col_mesh(ma)
            write(12,'(a)') trim(adjustl(snum))
            write(12,'(a)') trim(adjustl(scol))
            write(12,'(a)') 'gmt psxy _mesh${num}.txt -R -J -W0.5,#dddddd -G$col -K -O >> $fig'
        end do
        write(12,'(a)') 'gmt psxy _dis.txt -R -J -W0.5,#aaaaaa -K -O >> $fig'
        write(12,'(a)') 'echo 0 0 | gmt psxy -R -J -Sp -O >> $fig'
        write(12,'()')
!        write(12,'(a)') 'rm _*'
        close(12)
    end subroutine MAKE_BAT


    subroutine PST_PL(sigx,sigy,tauxy,ps1,ps2,ang)
        !Principal stress
        real(4),intent(in)::sigx,sigy,tauxy
        real(4),intent(out)::ps1,ps2,ang

        ps1=0.5*(sigx+sigy)+sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        ps2=0.5*(sigx+sigy)-sqrt(0.25*(sigx-sigy)*(sigx-sigy)+tauxy*tauxy)
        if(sigx==sigy)then
            if(tauxy>0.0)ang=45.0
            if(tauxy<0.0)ang=135.0
            if(tauxy==0.0)ang=0.0
        else
            ang=0.5*atan(2.0*tauxy/(sigx-sigy))
            ang=180.0/pi*ang
            if(sigx>sigy.and.tauxy>=0.0)ang=ang
            if(sigx>sigy.and.tauxy<0.0)ang=ang+180.0
            if(sigx<sigy)ang=ang+90.0
        end if
    end subroutine PST_PL


    real(4) function AREA(x1,y1,x2,y2,x3,y3)
        real(4),intent(in)::x1,y1,x2,y2,x3,y3
        AREA=0.5*((x3-x2)*y1+(x1-x3)*y2+(x2-x1)*y3)
    end function AREA

end program f90_gmt_ini4_vec
