program poscp
c------------------------------
c ce prog. permet de calculer la position (notée d)
c du centre de poussée par rapport au foyer 
c sur la corde du profil
c------------------------------
c.. déclaration des variables
	integer nb,ok
	character*72 nom,noms
	real*8 a,i,cz,cx,poub,cm,pi,pente,cz1,cz2,cm1,cm2,a1,a2,x,z,d
	parameter (pi=3.1415926)
	ok = 0
c.. ouverture des fichiers de données Xfoil (11) et de sortie (12)
	print*,'-- nom du fichier de données Xfoil contenant'
	print*,'   i Cz Cx ... Cm ...  ?'
	read*, nom
	open(unit=11,file=nom)
	do i=1,12
	  read (11,*)
	end do
	print*,'-- nombre d angles à gérer ?'
	read*, nb
	print*,'-- nom du fichier de sortie ?'
	read*, noms
	open(unit=12,file=noms)
c.. première boucle sur le fichier d'entrée Xfoil permettent de
c.. rechercher le coefficient de moment pour une portance nulle
c.. Cm0 = Cm(Cz=0)
	do j=1, nb
	  read (11,*) a,cz,poub,poub,cm
c	  print*, a,cz,poub,poub,cm
	  if (j .eq. 1) then
	    cz1 = cz
	    cm1 = cm
	    a1 = a
	  else
	    if ((cz*cz1) .lt. 0) then
	      print*,'** changement de signe de Cz entre',a1, 'et ',a
	      cz2 = cz
	      cm2 = cm
	      ok = 1
	      goto 50
	    else
	      cz1 = cz
  	      cm1 = cm
	      a1 = a  
	    endif
	  endif
	end do
50	continue
	if (ok .eq. 0) then
	  print*,'%% Il n y a pas de changement de signe sur Cz',
     &' dans votre fichier de données'
     	  stop  
	endif
c	print*,cm2,cm1,cz2,cz1
	pente = (cm2-cm1)/(cz2-cz1)
	cm0 = ( cm1 + cm2 - pente*(cz1+cz2) ) / 2.0
	print*,''
	print*,'** valeur de Cm0 = ',cm0
	print*,''
	rewind(11)
c.. nouvelle boucle sur le fichier d'entrée pour créer les courbes
c.. donnant les positions exactes et approchées du centre de poussée  
	do i=1,12
	  read (11,*)
	end do
	write (12,99)
	do j=1, nb
	  read (11,*) a,cz,cx,poub,cm
c	  print*, a,cz,cx,poub,cm
	  i = a * pi / 180.0
	  x = - cm / (cz+cx*tan(i))
	  z = - x * tan(i)
c.. distance par rapport au foyer
	  d = sqrt ( x**2 + z**2 ) * x / abs ( x )
c	  print*,x,z,d
c.. distance par rapport au bord d'attaque
	  d = d + 0.250
	  write (12,100) a,d,(0.250-cm0/cz),cx,cz,cz/cx,cm
c	  print 100 a,d,(0.250-cm0/cz),cx,cz,cz/cx,cm
	end do
	print *,'** fichier de sortie normalement écrit :',noms
c........1.........2.........3.........4.........5.........6.........7..
99      format('# i(°)        d_exact     d_approchee  Cx          Cz  
     &      Cz/Cx          Cm')
100	format(f9.4,3x,f9.4,3x,f9.4,3x,f9.4,3x,f9.4,3x,f9.4,3x,f9.4)
	end