! Ferma-cablu
! Cu Ecuatia neliniara a sagetii: 1; 2 3.
!  Fisierul de date:
!  l,f1,f2, n    : semi-deschidere, sageti; nr. de panouri
!  A1,A2,A3   : arii - cablu sup., cablu inf., montanti
!  H          : forta de pretensionare - c. inf.
!  E,sigma    : modul de elast.; tensiune limita
! iEq: indicele Ecuatiiei sagetii: 1; 2; >=3.
!
program NELINIAR
real l,l2, ln
character(260) fname, dummy*80
logical alta_incarcare
real(4) XX(2), F(2), JF(2,2)

common B1, B2, B3, C1, C2, c0, r21
common ieq

EXTERNAL fcn, jfcn
EXTERNAL sub1, sub2

character(1) char1

                  
1 print '(a, \)',' Numele fisierului de date: '
read '(a)', fname
open(1, file=fname, status='old')

90 read(1,'(a)') dummy		! Antet
if(dummy(1:1) =='/' .and. dummy(2:2) =='/') then
goto 90
else
backspace 1
endif

read(1,*) l, f10, f20, NP
read(1,*) A1, A2, A3
read(1,*) H20
read(1,*) E, sigma

read(1,*) ieq


l2 = l*l			!a = 2*l/n
! Pretensionare:
H10 = H20*f20/f10
p0  = H20*f20*2/l2
write(1,*)

! incarcarea initiala in montanti
print*,'  Incarcare echivalenta din pretensionare:', p0
write(1,*) '  Incarcare echivalenta din pretensionare:', p0

! Ecuatia:
write(1, '(/, a, I1,/)') 'Ecuatia sagetii - ', ieq


! Constante:
EA1 =E*A1; 
EA2 =E*A2;
r21 =A2/A1

x10 =f10/l; 
x20 =f20/l;
c0 =x10+x20

B10 =p0*l/EA1;
B20 =p0*l/EA2;

C1 =B10/x10 -phi(x10)
C2 =B20/x20 -phi(x20)


! Incarcarea de slabire montanti: dp=p0
SELECT CASE(ieq)

CASE(1)
! Solutii explicite
x2sm =0.5*sqrt( exp(-C2)-3 )
x1sm =c0 -x2sm

CASE(2:)
! Ecuatie: phi(x2) +C2=0
eps =1e-6; lnit =10; kod =0;

call Newton(SUB2,x10,eps,lnit,rad,kod)

x2sm =rad
x1sm =c0-x2sm

END SELECT


B1sm =x1sm*(phi(x1sm) +C1)
q_sl_m =B1sm*EA1/l

print*,' *Incarcare de slabire montanti          :', q_sl_m
write(1,*) ' *Incarcare de slabire montanti          :', q_sl_m

! incarcarea de slabire cablu superior: =OK
!q_sl_cs = -p
SELECT CASE(ieq)

CASE(1)
x1scs =0.5*sqrt( exp(-C1)-3 )
x2scs =c0 -x1scs


CASE(2:)
! Ecuatie: phi(x1) +C1=0
eps =1e-6; lnit =10; kod =0;

call Newton(SUB1, x10, eps, lnit, rad, kod)

x1scs =rad
x2scs =c0-x1scs

END SELECT


B2scs =x2scs*(phi(x2scs) +C2)
p_sl_cs =EA2*B2scs/l
q_sl_cs = -p_sl_cs

print*,' *Incarcare de slabire cablu superior    :', q_sl_cs
write(1,*) ' *Incarcare de slabire cablu superior    :', q_sl_cs
      
icaz=0
print'(a)',' Introdu incarcarea: q_min, q_max, q_pas (q, 0, 0 = pt. o singura incarcare):'
read*, q_min, q_max, q_pas
if(q_max ==0 .and. q_pas ==0) then
q_max =q_min
q_pas =1
endif


!  Write results header:
write(1,10)' Caz','Incarcare','Sageata', '"df"',   'H1',       'H2',      'T-montanti'
10 format(// t3,a, t12,a,      t26,a,       t40,a,    t52,a,      t64,a,     t76,a, /)
11 format(' ',I5,  t11,f10.3, t24,f10.6,   t36,f10.6, t49,f11.4,  t62,f10.4, t73, f10.4)

n =2
lnit0 =10
xtol =1e-6
xx(1) =x10; xx(2) =x20;
a =2*l/NP

100 DO q=q_min, q_max, q_pas
icaz = icaz+1
B3 =q*l/EA1
lnit =lnit0; kod =0;

call Newton_Sys(fcn, jfcn, n, xx, lnit, xtol, f, kod)


x1 =xx(1); x2 =xx(2)

B2 =x2*(phi(x2)+C2)

f1 =x1*l; f2 =x2*l
p =EA2*B2/l

H1 =(q+p)*l2/(2*f1)
H2 =p*l2/(2*f2)
!print*, ' p=',p
!write(1,*) ' p=',p
Tm =p*a


!  Write Results:
write(1,11) icaz, q, f1, f1-f10, H1, H2, Tm

if(H1/A1 >=sigma) then
	print*,' * Cablu superior rupt'
	write(1,*) ' * Cablu superior rupt: sigma =', H1/A1
endif

if(H1 <=0) then
	print*,' * Cablu superior slabit'
    write(1,*) ' * Cablu superior slabit'
endif

if(H2/A2 >=sigma) then
	print*,' * Cablu inf. rupt'
    write(1,*) ' * Cablu inferior rupt: sigma =', H2/A2
endif

if(H2 <=0) then
	print*,' * Cablu inf. slabit'
    write(1,*) ' * Cablu inferior slabit'
endif

if(p <=0.) then
	print*,' * Montanti slabiti'
    write(1,*) ' * Montanti slabiti'
elseif(Tm/A3 >=sigma)then
	print*,' * Montanti rupti'
    write(1,*) ' * Montanti rupti: sigma =', Tm/A3
endif
write(1,*)
ENDDO

print '(/, a, \)', ' Alta incarcare ? (Y, y / char): '
read(*, '(a1)', err =500) char1

alta_incarcare =char1 =='Y' .or. char1 =='y'

if(alta_incarcare) then
! Read "q":
print '(a)',' Introdu: q_min, q_max, q_pas (q, 0, 0 = pt. o singura incarcare): '

read *, q_min, q_max, q_pas
	if(q_max ==0 .and. q_pas ==0) then
	q_max =q_min
	q_pas =1
	endif
goto 100
endif

500   close (1)

print '(a, a, /)',' * Fisier de rezultate :', trim(fname)

END

