program bdforce4
c------------------------------
c ce prog. permet de tracer les bd d'évolution de force sur un profil
c à un nombre de Reynolds et cela à partir des fichiers
c issus de Xfoil comportant une série d'angle.
c en mettant des valeurs de force
c------------------------------
	integer nb,j,jj,pass
	character*72 nom
	character*6 justif
	real a(60),cz(60),cx(60),poub,cm(60),i,x,z,d,Re,c1,c2,
     &       dmin,dmax,ampv,pi,nu,rho,co,vi,coef
        parameter (pi=3.1415926 , nu=0.150e-4 , rho=1.20 )
	justif=" left "
	pass=1
c.. ouverture du fichier de données Xfoil
	print*,'-- nombre d angles a gerer ?'
	read*, nb
	print*,'-- nombre Re ?'
	read*, Re
	print*,'-- corde du profil ?'
	read*, co
	vi = nu*Re/co
	print*,'** Vitesse ',vi,' m/s'
	coef = 0.50*rho*co*vi**2
	  print*,'-- nom du fichier de donnees Xfoil contenant'
	  print*,'   i Cz Cx ... Cm ...  ?'
	  read*, nom
	  open(unit=11,file=nom)
	  do jj=1,12
	    read (11,*)
	  end do
	  do jj=1,nb
	    read (11,*) a(jj),cz(jj),cx(jj),poub,cm(jj)
	    print*,jj, a(jj),cz(jj),cx(jj),poub,cm(jj)
	  end do
c.. fin de lecture des informations Xfoil
c
	open (unit=13, file='arrow')
	dmin = -0.5
	dmax = 1.5
	ampv = (dmax-dmin)/10.0
	ampf = (dmax-dmin)/1.0
        write (13,112) Re,co,vi
112     format('set title "Re=',f9.0,' ( corde= ',f5.3,' m _',
     &         ' vitesse= ',f4.1,' m/s )"',/,
     &         'set noborder',/,
     &         'set noxtics',/,
     &         'set noytics',/,
     &         'set xzeroaxis',/,
     &         'set nokey',/,
     &         'set size 0.721,1.0',/,
c     &         'set term jpeg',/,
     &         'set term png',/,
     &         'set data style lines',/,
     &	       'set xrange [-0.5:1.5]',/,
     &         'set yrange [-0.5:1.5]')
c.. boucle pour représenter les vecteurs forces et vitesse sur chaque angle pour les 2 profils d'aile
	do jj=1,nb
	   print*,jj
c.. dessin remis vierge
	    write (13,120) 	
	    i = a(jj) * pi / 180.0
c valeur de l'angle
	    write (13,99) a(jj)
c vecteurs vitesse
	    write (13,103) 2,dmin,(-ampv),(dmin+ampv),((tan(i)-1.0)*ampv)
	    write (13,103) 3,dmin,0.0,(dmin+ampv),(tan(i)*ampv)
	    write (13,103) 4,dmin,ampv,(dmin+ampv),((tan(i)+1.0)*ampv)
c.. distance par rapport au bord d'attaque
	    d = 0.25 - cm(jj) / (cz(jj)*cos(i)+cx(jj)*sin(i))
c vecteur force
	      c1 = (cx(jj)*cos(i) - cz(jj)*sin(i))*ampf
	      if ((d+c1) .gt. 1.50 ) then
c.. vecteur sort du dessin
	        c1=1.50-d
	      endif
	      c2 = (cx(jj)*sin(i) + cz(jj)*cos(i))*ampf
	      if (c2 .gt. 1.50 ) then
c.. vecteur sort du dessin
	        c2 = 1.50
	      endif
	      if (c2 .lt. -0.50 ) then
c.. vecteur sort du dessin
	        c2 = -0.50
	      endif
	      if ( d .ge. -0.5 .and. d .le. 1.5) then
		write (13,100) 6,d,0,d+c1,c2,2
	      else
		write (13,98)   
	      endif
c.. valeur de la force
	      fo = sqrt( cx(jj)**2+cz(jj)**2 ) * coef
	      write (13,104) 2,fo,d+c1,c2,justif
c nom des fichiers jpeg 
	    write (13,101) (9+jj)
	    if (pass .eq. 1) then
	      write (13,131)
	      pass=2
	    else
	      write (13,132)
	    endif
	end do
	print*,'** creer dossier : mkdir bd'
	print*,'** copier fichier profil :'
	print*,'   * cp "votre profil" airfoil1.dat'
	print*,'   * les coordonnees du fichier profil commencent'
     &' sur la premiere ligne !'
	print*,'** lancer        : gnuplot<arrow'
120   format('set noarrow ',/,'set nolabel')
98    format('set label 2 "CP out"')
99    format('set label 1 "',f7.3,'°" at -0.4,1')
100   format('set arrow ',i1,' from ',f7.4,',',f7.4,' to ',f7.4,',',
     &       f7.4,' lt ',i1,' lw 2')
c101   format('set output "bd/',i2,'.jpg"')
101   format('set output "bd/',i2,'.png"')
131   format('plot "airfoil1.dat" lt 1')
132   format('replot')
103   format('set arrow ',i1,' from ',f7.4,',',f7.4,
     &         ' to ',f7.4,',',f7.4,' lt 3 lw 3')
104   format('set label ',i1,' "',f7.3,' N/m" at',f7.4,',',f7.4,a6)
	end