cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c LEPARACHUTE 0.12 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Vesrion 0.1 Date 2016-05-26 c Version 0.11 Date 2016-06-05 c Version 0.12 Date 2016-09-18 c Version 0.2 Date 2020-04-16 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Laboratoridenvol c http://www.laboratoridenvol.com cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c GNU General Public Licence 3.0 http://www.gnu.org cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc program leparachute ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 1. Variables declaration ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc character*50 wname, bname real*8 xwf,xkf real*8 a,b,c,ap,ar,lh,sep,calage real*8 x,y,xg,yg real*8 kparam, sfactor real*8 omega1,omega2,tetha,pi,asector,alpha,beta real*8 alpin,alpcir,alphag real*8 tlen(1000),dlen(1000),slen(1000) real*8 xp(0:1000),zp(0:1000),up(1000),vp(1000) real*8 upg(1000,1000),vpg(1000,1000) real*8 upgl(1000,1000),vpgl(1000,1000) real*8 x3(1000,1000),y3(1000,1000),z3(1000,1000) real*8 xzero,xrad,xradfi real*8 opx,opy,oqx,oqy real*8 xc1,yc1,xc2,yc2 real*8 xc1g(500,1000),yc1g(500,1000) real*8 xc2g(500,1000),yc2g(500,1000) real*8 bncxr,bncyr,bnczr,bncxl,bncyl,bnczl real*8 ancx(300,1000),ancy(300,1000),ancz(300,1000) real*8 xline(300,1000) real*8 areap1,areap2,areasail integer ngores,nmeri,npoin real*8 llarga,llargaa,llargax,llargam,llargad,llargas integer illarga real*8 upc,vpc,upcl,vpcl cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 2. Init cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Open units open(unit=20,file='parachute.dxf') open(unit=22,file='data.txt') open(unit=23,file='lepc-out.txt') open(unit=30,file='lines.txt') open(unit=25,file='lepc-3d.dxf') c Init dxf files call dxfinit(20) call dxfinit(25) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 3. Data reading ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc rewind (22) rewind (23) rewind (30) c Read data file do i=1,4 read (22,*) end do read (22,*) pcversion do i=1,4 read (22,*) end do read (22,*) bname ! Brand name read (22,*) read (22,*) wname ! Wing name read (22,*) read (22,*) xwf ! wing scale read (22,*) read (22,*) xkf ! drawing scale read (22,*) read (22,*) atp ! parachute type read (22,*) read (22,*) a ! a parameter read (22,*) read (22,*) b ! b parameter read (22,*) read (22,*) c ! c parameter read (22,*) read (22,*) ap ! ap parameter read (22,*) read (22,*) ar ! Apex radius read (22,*) read (22,*) lh ! Line height read (22,*) read (22,*) sep ! separation read (22,*) read (22,*) calage ! calage read (22,*) read (22,*) ngores ! gores read (22,*) read (22,*) sfactor ! general ovalization factor read (22,*) read (22,*) npoin ! npoints used in meridian lines c Change scale a=a*xwf b=b*xwf c=c*xwf ap=ap*xwf ar=ar*xwf lh=lh*xwf sep=sep*xwf calage=calage*xwf cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 4. Parameters cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc pi=4.*atan(1.0d0) c npoin=101 omega2=datan(c/ap) omega1=dacos(ap/a) kparam=(b*dsin(omega1)-c)/(ap*ap) tetha=pi-omega1 teinc=(pi-omega1)/float(npoin-1) nemeri=ngores asector=2.*pi/float(ngores) call evenodd(ngores,ieo) ! detect if even or odd areap1=(pi*(a+ap)**2.)/10000. ! area projected m2 areap2=areap1 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 5. Set meridian line ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,npoin if (tetha.le.(pi/2.)) then xg=a*dcos(tetha) zg=b*dsin(tetha) else xg=a*dcos(tetha) zg=b*dsin(tetha)-kparam*xg*xg end if xp(i)=xg+ap zp(i)=zg tetha=tetha-teinc end do ! i c Compute apex llarga=0.0d0 do i=1,npoin llarga=llarga+dsqrt((xp(i)-xp(i+1))**2.+((zp(i)-zp(i+1)))**2.) c write (*,*) i,llarga,xp(i),zp(i) c write (*,*) i,xp(i),xp(i+1),llarga if (xp(i).lt.ar.and.xp(i+1).ge.ar) then llargax=xp(i+1)-ar llargam=(zp(i+1)-zp(i))/(xp(i+1)-xp(i)) llargad=llargax*dsqrt(1+llargam*llargam) llargaa=llarga-llargad illarga=i llargas=dsqrt((xp(i)-xp(i+1))**2.+((zp(i)-zp(i+1)))**2.) c write (*,*) "> ",llargaa,llargad,illarga,llargas end if end do c Draw 2D dome do i=1,npoin-1 call line(xp(i),2000+zp(i),xp(i+1),2000+zp(i+1),4) call line(-xp(i),2000+zp(i),-xp(i+1),2000+zp(i+1),5) end do ! i cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 6. "t" lenght cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc tlen(1)=0.0d0 do i=2,npoin tlen(i)=tlen(i-1)+sqrt((xp(i)-xp(i-1))**2.+(zp(i)-zp(i-1)) + **2.) end do cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 7. "d" lenght and "s" lenght cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,npoin dlen(i)=2.*xp(i)*dsin(asector/2.) slen(i)=dlen(i)*(1.+sfactor/100.0d0) end do cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 8. Gore calculus cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calcule gore iteratively i=1 up(i)=0.0d0 vp(i)=0.0d0 do i=2,npoin up(i)=up(i-1)+0.5*(slen(i)-slen(i-1)) vp(i)=vp(i-1)-sqrt((tlen(i)-tlen(i-1))**2.-(slen(i)-slen(i-1)) + **2.) end do c Gore area area=0.0d0 do i=1,npoin-1 area=area+2.*(up(i)+up(i+1))*0.5*(vp(i)-vp(i+1)) end do areasail=area*float(ngores)/10000. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 9. Gore drawing (2D) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,npoin vp(i)=vp(i)-80.0d0 end do do k=1,ngores alphag=asector*float(k-1) do i=1,npoin upg(i,k)=up(i)*dcos(alphag)-vp(i)*dsin(alphag) vpg(i,k)=up(i)*dsin(alphag)+vp(i)*dcos(alphag) upgl(i,k)=-up(i)*dcos(alphag)-vp(i)*dsin(alphag) vpgl(i,k)=-up(i)*dsin(alphag)+vp(i)*dcos(alphag) end do ! i end do ! k do k=1,ngores do i=1,npoin-1 call line(upg(i,k),vpg(i,k),upg(i+1,k),vpg(i+1,k),1) call line(upgl(i,k),vpgl(i,k),upgl(i+1,k),vpgl(i+1,k),1) end do ! i i=npoin call line(upg(i,k),vpg(i,k),upgl(i,k),vpgl(i,k),1) end do ! k c Set cut for apex i=illarga+1 do k=1,ngores upc=upg(i,k)+(upg(i+1,k)-upg(i,k))*(llargas-llargad)/llargas vpc=vpg(i,k)+(vpg(i+1,k)-vpg(i,k))*(llargas-llargad)/llargas upcl=upgl(i,k)+(upgl(i+1,k)-upgl(i,k))*(llargas-llargad)/llargas vpcl=vpgl(i,k)+(vpgl(i+1,k)-vpgl(i,k))*(llargas-llargad)/llargas call line(upc,vpc,upcl,vpcl,1) c write (*,*) k,upc,vpc,upcl,vpcl end do ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 10. Gore drawing (3D) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c For each gore do k=1,ngores alpha=asector*float(k-1)+asector*0.5 do i=1,npoin x3(i,k)=xp(i)*dcos(alpha) y3(i,k)=xp(i)*dsin(alpha) z3(i,k)=zp(i) end do c Meridian line do i=1,npoin-1 call line3d(x3(i,k),y3(i,k),z3(i,k),x3(i+1,k),y3(i+1,k), + z3(i+1,k),3) end do end do ! end k c Draw straight segments (optional) i=npoin do k=1,ngores-1 call line3d(x3(i,k),y3(i,k),z3(i,k),x3(i,k+1),y3(i,k+1), + z3(i,k+1),3) end do k=ngores call line3d(x3(i,1),y3(i,1),z3(i,1),x3(i,k),y3(i,k), + z3(i,k),3) c Draw ovalized arc segments in 3D only for reference i=npoin xrad=0.5*dlen(i) c (IMPROVE) METHOD TO SOLVE USING BRUTE FORCE c Solve equation do j=1,1000 xzero=0.5*dlen(i)-xrad*dsin(0.5*slen(i)/xrad) if (abs(xzero).le.0.1) then xradfi=xrad end if xrad=xrad+2.*(ap+a)/1000 end do alpha=dlen(i)/(2.*xradfi) k=1 oqx=(x3(i,k)+x3(i,k+1))/2. oqy=(y3(i,k)+y3(i,k+1))/2. beta=abs(datan((x3(i,k+1)-x3(i,k))/(y3(i,k+1)-y3(i,k)))) opx=oqx-xradfi*dcos(alpha)*dcos(beta) opy=oqy-xradfi*dcos(alpha)*dsin(beta) c Recalcule xradfi xradfi=sqrt((x3(i,k)-opx)**2.+(y3(i,k)-opy)**2.) c Reacalcule alpha alpha=dasin(dlen(i)/(2.*xradfi)) c For each gore do k=1,ngores alphag=asector*float(k-1) c Arc alpin=alpha*2./30. alpcir=0. do j=1,30 xc1=(opx+xradfi*dcos(alpha+beta+alpcir)) yc1=(opy+xradfi*dsin(alpha+beta+alpcir)) xc2=(opx+xradfi*dcos(alpha+beta+alpcir-alpin)) yc2=(opy+xradfi*dsin(alpha+beta+alpcir-alpin)) xc1g(j,k)=xc1*dcos(alphag)-yc1*dsin(alphag) yc1g(j,k)=xc1*dsin(alphag)+yc1*dcos(alphag) xc2g(j,k)=xc2*dcos(alphag)-yc2*dsin(alphag) yc2g(j,k)=xc2*dsin(alphag)+yc2*dcos(alphag) call line3d(xc1g(j,k),yc1g(j,k),z3(i,k),xc2g(j,k),yc2g(j,k), + z3(i,k),3) alpcir=alpcir-alpin end do end do ! k ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 11. LINES ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Main anchors bncxr=calage bncyr=sep/2. bnczr=-lh bncxl=calage bncyl=-sep/2. bnczl=-lh c Anchor points i=npoin do k=1,ngores ancx(1,k)=x3(i,k) ancy(1,k)=y3(i,k) ancz(1,k)=z3(i,k) end do c Lines c Case even if (ieo.eq.2) then do k=1,ngores/2 call line3d(bncxr,bncyr,bnczr,ancx(1,k),ancy(1,k),ancz(1,k),9) xline(1,k)=sqrt((bncxr-ancx(1,k))**2.+(bncyr-ancy(1,k))**2.+ + (bnczr-ancz(1,k))**2.) end do do k=1+ngores/2,ngores call line3d(bncxl,bncyl,bnczl,ancx(1,k),ancy(1,k),ancz(1,k),9) xline(1,k)=sqrt((bncxl-ancx(1,k))**2.+(bncyl-ancy(1,k))**2.+ + (bnczl-ancz(1,k))**2.) end do end if c Caso odd if (ieo.eq.1) then do k=1,(ngores+1)/2 call line3d(bncxr,bncyr,bnczr,ancx(1,k),ancy(1,k),ancz(1,k),9) xline(1,k)=sqrt((bncxr-ancx(1,k))**2.+(bncyr-ancy(1,k))**2.+ + (bnczr-ancz(1,k))**2.) end do do k=(ngores+1)/2,ngores call line3d(bncxl,bncyl,bnczl,ancx(1,k),ancy(1,k),ancz(1,k),9) xline(1,k)=sqrt((bncxl-ancx(1,k))**2.+(bncyl-ancy(1,k))**2.+ + (bnczl-ancz(1,k))**2.) end do end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Write report cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Print in report file write (23,*) "PARACHUTE AND PARASAIL DESIGN PROGRAM" write (23,*) "Laboratori d'envol" write (23,'(A25,F6.2)') "Version: ",pcversion write (23,'(A25,A50)') "Brand: ", bname write (23,'(A25,A50)') "Model: ", wname write (23,'(A25,i4)') "N gores: ", ngores write (23,'(A25,F7.2)') "Main radius (m): ",(a+ap)/100. write (23,'(A25,F7.2)') "Dome height (m): ",lh/100. write (23,'(A25,F7.2)') "Meridian lenght (m): ",tlen(npoin)/100. write (23,'(A25,F7.2)') "Surface sail (m2): ",areasail write (23,'(A25,F7.2)') "Surface projection (m2): ",areap1 c Print console write (*,'(A25,i4)') "N gores: ",ngores write (*,'(A25,F7.2)') "Main radius (m): ",(a+ap)/100. write (*,'(A25,F7.2)') "Dome height (m): ",lh/100. write (*,'(A25,F7.2)') "Meridian lenght (m): ",tlen(npoin)/100. write (*,'(A25,F7.2)') "Surface sail (m2): ",areasail write (*,'(A25,F7.2)') "Surface projection (m2): ",areap1 write (*,*) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Write lines cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc write (23,*) write (23,'(A14)') "Lines (cm): " write (30,'(A14)') "Lines (cm): " do k=1,ngores write (23,'(A8,I3,F10.2)') "Line ",k,xline(1,k) write (30,'(A8,I3,F10.2)') "Line ",k,xline(1,k) end do write (23,'(A15,F10.2)') "Central Line ",lh+c write (30,'(A15,F10.2)') "Central Line ",lh+c write (23,'(A21,F10.2)') "Central apex lines ",llargaa write (30,'(A21,F10.2)') "Central apex lines ",llargaa ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c END MAIN PROGRAM ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (atp.eq.1) then write (*,*) "OK Parachute calculated" end if if (atp.eq.2) then write (*,*) "OK Parasail calculated" end if call dxfend(20) call dxfend(25) close(20) close(22) close(23) close(30) close(25) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINES ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE LINE 2D ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE line(p1x,p1y,p2x,p2y,linecolor) c line P1-P2 real*8 p1x,p1y,p2x,p2y write(20,'(A,/,I1,/,A)') "LINE",8,"default" write(20,'(I1,/,A)') 6,"CONTINUOUS" write(20,'(I2,/,F14.4,/,I2,/,F14.4)') 10,p1x,20,p1y write(20,'(I2,/,F14.4,/,I2,/,F14.4)') 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 real*8 p1x,p1y,p1z,p2x,p2y,p2z write(25,'(A,/,I1,/,A)') "LINE",8,"default" c write(25,'(I3,/,A)') 100,"AcDbLine" write(25,'(I1,/,A)') 6,"CONTINUOUS" write(25,'(I2,/,F14.3,/,I2,/,F14.3,/,I2,/,F14.3)') + 10,p1x,20,p1y,30,p1z write(25,'(I2,/,F14.3,/,I2,/,F14.3,/,I2,/,F14.3)') + 11,p2x,21,p2y,31,p2z write(25,'(I2,/,I2,/,I2,/,I2,/,I2)') 39,0,62,linecolor,0 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,-670.,20,-3630. write(nunit,'(I1,/,A)') 9,"$EXTMIN" write(nunit,'(I2,/,F12.3,/,I2,/,F12.3)') 10,7000.,20,120. 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 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c EVEN (parell chetnoye) ODD (senar nechetnoye) subroutine ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine evenodd(ngores,ieo) control=((ngores)/2.)-float(int((ngores)/2.)) if (control.ne.0) then ieo=1 else ieo=2 end if return end