c***************************************************************
c      ASN AIRFOIL TO SHARK NOSE AIRFOIL MODIFICATION
c      Convert any .dat profile into a shark-nose profile using 
c      modification functions with two generic parabolas of degree n
c      Pere Casellas
c      Laboratori d'envol  
c      https://www.laboratoridenvol.com
c      v0.1 20260802
c      GNU General Public License 3.0 (http://www.gnu.org)
c***************************************************************
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      VARIABLES DECLARATION
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       real x1,y1,z1,x2,y2,z2
       real p1x,p1y,p1z,p2x,p2y,p2z
       real xmod1,xmod2,xmod3
       real xdfle1,xdfle3,xn1,xn3
       real xxco(1000), xyco(1000)
       character(50) aname, gname

c       pi=4.*atan(1.)

       write (*,*) 
       write (*,*) "**************************************************"
       write (*,*) "PROGRAM airfoil shark nose modification"
       write (*,*) "verson 0.1 20260802"
       write (*,*)
       write (*,*) "Pere Casellas"
       write (*,*)
       write (*,*) "GNU General Public License 3.0 http://www.gnu.org"
       write (*,*) "**************************************************"
       write (*,*)
       write (*,*)
      

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Open files
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       open(unit=20, file='airfoil-sn.dxf')
       open(unit=21, file='asndata.txt')
       open(unit=22, IOSTAT=nios, file='airfoil.dat')
       open(unit=23,file='airfoil-sn.dat')
       
       call dxfinit(20)
       
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Reading data file asndata.txt
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
       rewind (21)
       rewind (22)
       rewind (23)

c      Read data

       read (21,*) xmod1,xdfle1,xn1
       read (21,*) xmod2
       read (21,*) xmod3,xdfle3,xn3
       read (21,*) xcale

       xmod1=xmod1*0.01
       xmod2=xmod2*0.01
       xmod3=xmod3*0.01
       xdfle1=xdfle1*0.01
       xdfle3=xdfle3*0.01
       imod1=0
       imod2=0
       imod3=0
       kplus=0

       k=0
       
c      Maximal allowed lines in file .dat = 10000  

c      Read airfoil name 
        
       read (22,*) aname
       write (*,*) "Airfoil name = ", aname

c      Read airfoil coordinates       
       
       do
       if (nios.ge.0) then           
       k=k+1            
       read (22,*,IOSTAT=nios) xxco(k), xyco(k)
       end if
       if (nios.lt.0) then
               exit
       end if
       end do
       kmax=k

       write (*,*)
       
       do k=1,kmax-1
       write (*,'(I4,2x,F9.6,4x,F9.6)') k,xxco(k), xyco(k)
       end do

cccccccccccccccccccccccccccccccccccccccccccccccc
c      Detect main points
cccccccccccccccccccccccccccccccccccccccccccccccc

       write (*,*) 
       write (*,*) "Main points:"
       write (*,*) 
       write (*,'(x,A3,2x,I4)') "kT=", kmax-1

       kplus=0 ! additional airfoil points

       do k=1,kmax-1

c      Nose point
       if (xxco(k).lt.0.00001.and.xyco(k)*xyco(k).lt.0.00001) then
       k0=k
       write (*,'(x,A3,2x,I4)') "k0=",k0
       end if

c      Fisrt modification point
       if (xxco(k).le.xmod1.and.xxco(k+1).gt.xmod1.
     + and.xyco(k).le.0.0) then
       k1=k
       write (*,'(x,A3,2x,I4,2x,F7.3,2x,F7.3,2x,F7.3)') 
     + "k1=",k1,xxco(k),xmod1,xxco(k+1)
c      Set interpolation point if necessary
       if ((xmod1-xxco(k)).gt.0.001) then
       kplus=kplus+1
       imod1=1
       xxco1=xmod1
       xyco1=xyco(k)+((xyco(k+1)-xyco(k))*(xmod1-xxco(k))/
     + (xxco(k+1)-xxco(k)))
       end if
       end if

c      Second modification point
       if (xxco(k).le.xmod2.and.xxco(k+1).gt.xmod2.
     + and.xyco(k).le.0.0) then
       k2=k
       write (*,'(x,A3,2x,I4,2x,F7.3,2x,F7.3,2x,F7.3)') 
     + "k2=",k2,xxco(k),xmod2,xxco(k+1)
c      Set interpolation point if necessary
       if ((xmod2-xxco(k)).gt.0.001) then
       kplus=kplus+1
       imod2=1
       xxco2=xmod2
       xyco2=xyco(k)+((xyco(k+1)-xyco(k))*(xmod2-xxco(k))/
     + (xxco(k+1)-xxco(k)))
       end if
       end if

c      Third modification point
       if (xxco(k).le.xmod3.and.xxco(k+1).gt.xmod3.
     + and.xyco(k).le.0.0) then
       k3=k
       write (*,'(x,A3,2x,I4,2x,F7.3,2x,F7.3,2x,F7.3)') 
     + "k3=",k3,xxco(k),xmod3,xxco(k+1)
c      Set interpolation point if necessary
       if ((xmod3-xxco(k)).gt.0.001) then
       kplus=kplus+1
       imod3=1
       xxco3=xmod3
       xyco3=xyco(k)+((xyco(k+1)-xyco(k))*(xmod3-xxco(k))/
     + (xxco(k+1)-xxco(k)))
       end if
       end if

       end do

       write (*,*)
       write (*,*) "Points modification by interpolation:"
       write (*,*) "Imod1: ", imod1
       write (*,*) "Imod2: ", imod2
       write (*,*) "Imod3: ", imod3
       write (*,*) "Additional points: ",kplus

cccccccccccccccccccccccccccccccccccccccccccccccc
c      Shark nose modifications
cccccccccccccccccccccccccccccccccccccccccccccccc

       do k=1,kmax

c      First modification
c      General parabola beetween point k0 and k1 
       if (xxco(k).ge.0.0.and.xxco(k).le.xmod1.and.xyco(k).le.0.0) then
       xyco(k)=xyco(k)+(xdfle1/(xmod1**xn1))*xxco(k)**xn1
       end if

c      Second modification beetween point k1 and k3 
       if (xxco(k).gt.xmod1.and.xxco(k).le.xmod3.
     + and.xyco(k).le.0.0) then
       xyco(k)=xyco(k)+(xdfle3/((xmod3-xmod1)**xn3))*
     + (xmod3-xxco(k))**xn3
       end if

       end do

c      Additional interpolated points if necessary
       if (imod1.eq.1) then
       xyco1=xyco1+(xdfle1/(xxco1**xn1))*xxco1**xn1
       end if
       if (imod2.eq.1) then
       xyco2=xyco2+(xdfle3/((xmod3-xmod1)**xn3))*
     + (xmod3-xxco2)**xn3
       end if
       if (imod3.eq.1) then
       xyco3=xyco3+(xdfle3/((xmod3-xmod1)**xn3))*
     + (xmod3-xxco3)**xn3
       end if


ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Write new file airfoil-sn.dat
ccccccccccccccccccccccccccccccccccccccccccccccccccccc

       gname=trim(aname)//"-SN"

       write (23,'(A50)') gname

       do k=1,k0
       write (23,'(F9.6,4x,F9.6)') xxco(k), xyco(k)*xcale
       end do

       do k=k0+1,k1
       write (23,'(F9.6,4x,F9.6)') xxco(k), xyco(k)*xcale
       end do

       if (imod1.eq.1) then
c       write (23,*)
       write (23,'(F9.6,4x,F9.6)') xxco1, xyco1*xcale
c       write (23,*)
       end if

       do k=k1+1,k2
       write (23,'(F9.6,4x,F9.6)') xxco(k), xyco(k)*xcale
       end do

       if (imod2.eq.1) then
c       write (23,*)
       write (23,'(F9.6,4x,F9.6)') xxco2, xyco2*xcale
c       write (23,*)
       end if

       do k=k2+1,k3
       write (23,'(F9.6,4x,F9.6)') xxco(k), xyco(k)*xcale
       end do

       if (imod3.eq.1) then
c       write (23,*)
       write (23,'(F9.6,4x,F9.6)') xxco3, xyco3*xcale
c       write (23,*)
       end if

       do k=k3+1,kmax-1
       write (23,'(F9.6,4x,F9.6)') xxco(k), xyco(k)*xcale
       end do


ccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Draw airfoil-sn.dxf
ccccccccccccccccccccccccccccccccccccccccccccccccccccc

c      Upper part (green line, color 3)
       do k=1,k0-1
       call line(100.*xxco(k), -100.*xyco(k)*xcale, 
     + 100.*xxco(k+1), -100.*xyco(k+1)*xcale, 3)
       end do

c      First modification (red line, color 1)
       do k=k0,k1-1
       call line(100.*xxco(k), -100.*xyco(k)*xcale, 
     + 100.*xxco(k+1), -100.*xyco(k+1)*xcale, 1)
       if (imod1.eq.0) then
       call line(100.*xxco(k1),-100.*xyco(k1)*xcale,100.*xxco(k1),0.0,7)
       end if
       end do

c      If interpolated point for mod1
       if (imod1.eq.1) then
       call line(100.*xxco(k1), -100.*xyco(k1)*xcale, 
     + 100.*xxco1, -100.*xyco1*xcale, 2)
       call line(100.*xxco1, -100.*xyco1*xcale, 
     + 100.*xxco(k1+1), -100.*xyco(k1+1)*xcale, 1)
       call line(100.*xxco1,-100.*xyco1*xcale,100.*xxco1,0.0,7)
       k1=k1+1
       end if
      
c      Inlet (cyan line, color 4)
       do k=k1,k2-1
       call line(100.*xxco(k), -100.*xyco(k)*xcale, 
     + 100.*xxco(k+1), -100.*xyco(k+1)*xcale, 4)
       if (imod2.eq.0) then
       call line(100.*xxco(k2),-100.*xyco(k2)*xcale,100.*xxco(k2),0.0,7)
       end if
       end do

c      If interpolated point for inlet
       if (imod2.eq.1) then
       call line(100.*xxco(k2), -100.*xyco(k2)*xcale, 
     + 100.*xxco2, -100.*xyco2*xcale, 2)
       call line(100.*xxco2, -100.*xyco2*xcale, 
     + 100.*xxco(k2+1), -100.*xyco(k2+1)*xcale, 1)
       call line(100.*xxco2,-100.*xyco2*xcale,100.*xxco2,0.0,7)
       k2=k2+1
       end if

c      Second modification (blue line, color 5)
       do k=k2,k3-1
       call line(100.*xxco(k), -100.*xyco(k)*xcale, 
     + 100.*xxco(k+1), -100.*xyco(k+1)*xcale, 5)
       if (imod3.eq.0) then
       call line(100.*xxco(k3),-100.*xyco(k3)*xcale,100.*xxco(k3),0.0,7)
       end if
       end do

c      If interpolated point for mod3
       if (imod3.eq.1) then
       call line(100.*xxco(k3), -100.*xyco(k3)*xcale, 
     + 100.*xxco3, -100.*xyco3*xcale, 2)
       call line(100.*xxco3, -100.*xyco3*xcale, 
     + 100.*xxco(k3+1), -100.*xyco(k3+1)*xcale, 1)
       call line(100.*xxco3,-100.*xyco3*xcale,100.*xxco3,0.0,7)
       k3=k3+1
       end if

c      Not modified airfoil (red line, color 3)
       do k=k3, kmax-1
       call line(100.*xxco(k), -100.*xyco(k)*xcale, 
     + 100.*xxco(k+1), -100.*xyco(k+1)*xcale, 3)
       end do

c      Reset k values
       if (imod1.eq.1) then
       k1=k1-1
       end if
       if (imod2.eq.1) then
       k2=k2-1
       end if
       if (imod3.eq.1) then
       k3=k3-1
       end if

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      Final del programa principal
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c       call dxfend(25)
       call dxfend(20)
       
       write (*,*)
       write (*,*) 
       write (*,*) "OK, airfoil shark nose modified "
       write (*,*) "Total point in new shark nose airfoil: ",
     + kmax-1+kplus
       write (*,*) "Open files: airfoil-sn.dat and airfoil-sn.dxf"
       write (*,*)
       close (20)
       close (21)
       close (22)
       close (23)
       close (25)


       end 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      SUBROUTINE LINE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       SUBROUTINE line(p1x,p1y,p2x,p2y,linecolor)
c      line P1-P2

       real x1,x2,y1,y2,z1,z2

       write(20,'(A,/,I1,/,A)') "LINE",8,"default"
       write(20,'(I1,/,A)') 6,"CONTINUOUS"
       write(20,'(I2,/,F12.2,/,I2,/,F12.2)') 10,p1x,20,-p1y
       write(20,'(I2,/,F12.2,/,I2,/,F12.2)') 11,p2x,21,-p2y
       write(20,'(I2,/,I2,/,I2,/,I2,/,I2)') 39,0,62,linecolor,0
       return
       end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc        
c     SUBROUTINE LINE 3D
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       SUBROUTINE line3d(p1x,p1y,p1z,p2x,p2y,p2z,linecolor)
c      line P1-P2
       write(25,'(A,/,I1,/,A)') "LINE",8,"default"
       write(25,'(I1,/,A)') 6,"CONTINUOUS"
       write(25,'(I2,/,F12.2,/,I2,/,F12.2,/,I2,/,F12.2)') 
     + 10,p1x,20,p1y,30,p1z
       write(25,'(I2,/,F12.2,/,I2,/,F12.2,/,I2,/,F12.2)') 
     + 11,p2x,21,p2y,31,p2z
       write(25,'(I2,/,I2,/,I2,/,I2,/,I2)') 39,0,62,linecolor,0
       return
       end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      POLYLINE 2D
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       SUBROUTINE poly2d(plx,ply,nvertex,linecolor)

       real plx(300),ply(300),plz(300)

       write (20,'(A,/,I1,/,I1,/,I2)') "POLYLINE",8,0,62
       write (20,'(I3,/,I2,/,I1)') linecolor,66,1
       write (20,'(I2,/,F3.1,/,I2,/,F3.1,/,I2,/,F3.1,/,I1)') 
     + 10,0.0,20,0.0,30,0.0,0
       
       do k=1,nvertex

       write (20,'(A,/,I1,/,I1,/,I2)') "VERTEX",8,0,62
       write (20,'(I3,/,I2,/,I1)') linecolor,66,1 
       write (20,'(I2,/,F9.3,/,I2,/,F9.3,/,I2,/,F9.3,/,I1)') 
     + 10,plx(k),20,ply(k),30,0.0,0

       end do
       
       write (20,'(A,/,I1,/,I1,/,I2)') "SEQEND",8,0,62
       write (20,'(I3,/,I1)') linecolor,0 

       write (*,'(A,/,I1,/,I1,/,I2)') "POLYLINE",8,0,62
       write (*,'(I3,/,I2,/,I1)') linecolor,66,1
       write (*,'(I2,/,F3.1,/,I2,/,F3.1,/,I2,/,F3.1,/,I1)') 
     + 10,0.0,20,0.0,30,0.0,0
       
       do k=1,nvertex

       write (*,'(A,/,I1,/,I1,/,I2)') "VERTEX",8,0,62
       write (*,'(I3,/,I2,/,I1)') linecolor,66,1 
       write (*,'(I2,/,F9.3,/,I2,/,F9.3,/,I2,/,F9.3,/,I1)') 
     + 10,plx(k),20,ply(k),30,0.0,0

       end do

       write (*,'(A,/,I1,/,I1,/,I2)') "SEQEND",8,0,62
       write (*,'(I3,/,I1)') linecolor,0 

       return

       end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      ELLIPSE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       SUBROUTINE ellipse(u0,v0,a,b,tet0,linecolor)

       real xe(300),ye(300)

       pi=4.*atan(1.)

       do ll=1,40

       tet=2.*pi*((float(ll)-1.)/39.)

c       write (*,*) ll,float(ll),tet," ",pi,"---"

       xe(ll)=u0+a*cos(tet)*cos(tet0)-b*sin(tet)*sin(tet0)
       ye(ll)=v0+a*cos(tet)*sin(tet0)+b*sin(tet)*cos(tet0)

       end do

       do ll=1,39

       p1x=xe(ll)
       p2x=xe(ll+1)
       p1y=ye(ll)
       p2y=ye(ll+1)

       call line(p1x,p1y,p2x,p2y,linecolor)

c       write (*,*) ll,tet*180./pi,xe(ll),ye(ll)

       end do

       return

       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c DXF init
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       
       SUBROUTINE dxfinit(nunit)
       
       write(nunit,'(I1,/,A,/,I1)') 0,"SECTION",2
       write(nunit,'(A)') "HEADER"
       write(nunit,'(I1,/,A)') 9,"$EXTMAX"
       write(nunit,'(I2,/,F12.3,/,I2,/,F12.3)') 10,-900.,20,90.
       write(nunit,'(I1,/,A)') 9,"$EXTMIN"
       write(nunit,'(I2,/,F12.3,/,I2,/,F12.3)') 10,5000.,20,-3000.
       write(nunit,'(I1,/,A,/,I1)') 0,"ENDSEC",0
       write(nunit,'(A,/,I1)') "SECTION",2
       write(nunit,'(A,/,I1)') "ENTITIES",0

       return
       end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c DXF end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE dxfend(nunit)

      write(nunit,'(A,/,I1,/,A)') "ENDSEC",0,"EOF"
      return
      end


      
