10 ! Program pentru fitarea unui spectru Mossbauer prin metoda
20 !   celor mai mici patrate  folosind metoda Gauss-Newton
30 !             cu algoritmul Levenberg-Marquart
40 !                     D. CIURCHEA 1993
50 !**************************************************************
60 !*             P R O G R A M U L  P R I N C I P A L                       *
70 !**************************************************************
80 RANDOMIZE(.75)
90 CALL CLEAR
100 !CREEN 2
110 !
120 !-------------    Specificarea variabilelor    ---------------
130 !
140 AA=1!Incepe evaluarea preciziei masinii
150 FOR I=1 TO 100
160 REPS=2^(-I)
170 IF AA=AA+REPS THEN 190
180 NEXT I
190 REPS=REPS*2 :: PRINT REPS :: ! Lungimea mantisei masinii.
200 FOR I=1 TO 70 :: A$=A$&" " :: NEXT I
210 NPCT=121! Nr. punctelor experimentale.
220 NPAR=5! nr. parametrilor din modelul teoretic.
230 DIM D(121)! Vector care contin spectrul experimental.
240 DIM DC(121)! Vector care contine spectrul calculat.
250 DIM W(5),V(5),VT(5)! Vect. pt. param. modelului.
260 DIM A(5,5),X(5),B(5),G(5)! Rez. sist. liniar.
270 DIM AINV(5,5)! Calculul matricii inverse.
280 NSIG=5! Nr. cifrelor semn. in raportarea rez. final.
290 PREC=10^(-NSIG)
300 MAXFN=5000! Nr. max. de evaluari ale func. permis de util.
310 PRINT "FITAREA UNUI SPECTRU MOSSBAUER SIMULAT"
320 PRINT "Programul este in executie"
330 DATA 575,8.4,85,39,60
340 DATA 480,14,70,52,73
350 SZ=0
360 FOR I=1 TO NPAR!  Citeste param. adevarati ai spectrului
370 READ V(I)!  experimental ai nitroprusiatului de Na.
380 NEXT I
390 GG=V(2)^2/4! Patratul semilargimii  la semiinaltime
400 MI=1E+20
410 MA=-1/MI
420 FOR I=1 TO NPCT
430 ZG=.2*(RND-.5)! Simularea erorilor experimentale.
440 BM=I-V(4)! Pozitia picului de la viteze negative,
450 BP=I-V(5)*2+V(4)!respectiv pt.viteze pozitive
460 !                          la centrul spectrului dat de V(5).
470 AA=V(2)*V(3)
480 !              AA= Aria picului,V(3) * semilargimea V(2) * 2.
490 D(I)=V(1)-AA/(BM*BM+GG)-AA/(BP*BP+GG)
500 !                 Vectorul D() contine intensitatea adevarata.
510 ZG=ZG*SQR(D(I))
520 !       Erorile "experimentale" urmeaza o distributie Poisson.
530 SZ=SZ+ZG*ZG! Suma erorilor introduse in spectru.
540 D(I)=D(I)+ZG! Intens. "masurata", cu er. "exper.".
550 IF MA<D(I)THEN MA=D(I)! Scalarea vectorului D().
560 IF MI>D(I)THEN MI=D(I)
570 NEXT I
580 SCL=(MA-MI)! Factorul de scala.
590 ALM=0!10!100 ! Parametrul Levenberg-Marquart.
600 FOR I=1 TO NPAR! Valorile initiale ale parametrilor.
610 READ V(I)
620 NEXT I
630 EPS=SQR(REPS)!*32! Cst. pt. evaluarea derivatei numerice.
640 !
650 !------------------       INITIALIZARE       -----------------
660 !
670 FOR I=1 TO NPAR! Initializarea vectorului de lucru.
680 W(I)=V(I)
690 NEXT I
700 GOSUB 2480!Valorarea functiei cu parametrii initiali.
710 F0=F
720 !
730 !--------------   INCEPUTUL PROCESULUI ITERATIV  -------------
740 !
750 !
760 !=======  CONSTRUIREA SISTEMULUI LINIAR DE ECUATII  ==========
770 !
780 FOR I=1 TO NPAR! Init. mat. Hessian si a term. liberi.
790 FOR J=I TO NPAR
800 A(I,J)=0
810 A(J,I)=0
820 NEXT J
830 W(I)=V(I)! Actualizarea vectorului de lucru.
840 B(I)=0
850 NEXT I
860 !
870 !------ Calculul derivatelor numerice  in fiecare punct ------
880 !
890 GAM=V(2)*V(2)/4
900 FOR IP=1 TO NPCT
910 F1=V(1)-V(2)*V(3)/((IP-V(4))^2+GAM)
920 F1=F1-V(3)*V(2)/((IP-V(5)*2+V(4))^2+GAM)
930 ! F1 este o referinta pentru calculul derivatelor numerice
940 DIF=D(IP)-F1
950 FOR I=1 TO NPAR!Calc. comp.gradientului in punctul IP.
960 AA=ABS(V(I)) :: IF AA<.1 THEN AA=.1! Scala pas.
970 PAS=AA*EPS :: W(I)=V(I)+PAS! Pas  deriv.num.
980 GAMW=W(2)*W(2)/4
990 F2=W(1)-W(2)*W(3)/((IP-W(4))^2+GAMW)
1000 F2=F2-W(3)*W(2)/((IP-W(5)*2+W(4))^2+GAMW)
1010 G(I)=(F2-F1)/PAS! Componenta gradientului.
1020 W(I)=V(I)! Refacerea vectorului de lucru.
1030 NEXT I
1040 !
1050 !--- Calculul matricii Hessian, A si al term. liberi, B  --
1060 !
1070 FOR I=1 TO NPAR
1080 FOR J=I+1 TO NPAR
1090 A(I,J)=A(I,J)+G(I)*G(J)! Mat. Hessian.
1100 A(J,I)=A(I,J)
1110 NEXT J
1120 A(I,I)=A(I,I)+G(I)*G(I)
1130 B(I)=B(I)+DIF*G(I)! Termenii liberi.
1140 NEXT I
1150 NEXT IP
1160 IFN=IFN+(NPAR*NPAR+NPAR)/2
1170 GOSUB 2640! Rezolvarea sist.liniar printr-o metoda adecvata.
1180 !
1190 !===================  AJUSTAREA PASULUI  ====================
1200 !
1210 IF INEF>6 OR IFN>MAXFN THEN 2060
1220 ALFA=1! Incercarea cu un pas cit mai mare
1230 GOSUB 3140
1240 F2=F
1250 IF F2>=F0 THEN 1290
1260 AL=ALFA
1270 F0=F2
1280 GOTO 1540
1290 ALFA=.5! Incercarea cu jumatatea pasului maxim
1300 GOSUB 3140
1310 F1=F
1320 IF F1>=F0 THEN 1360
1330 AL=ALFA
1340 F0=F1
1350 GOTO 1540
1360 ALFA=(3*F0-4*F1+F2)/4/(F0-F1-F1+F2)
1370 !       Ajust. pas. cu  o parabola prin abscisele  0, 0.5, 1.
1380 GOSUB 3140
1390 IF F>=F0 THEN 1440
1400 AL=ALFA
1410 F0=F
1420 GOTO 1540
1430 !
1440 !--------   Decizie pentru cautarea ineficienta   -----------
1450 !
1460 !IF F0 >= F1 AND F0 >= F2 THEN 1270
1470 ALM=ALM*1.01! Se amplifica param. L-M.
1480 INEF=INEF+1
1490 GOSUB 2410! Schimba semnul si marimea corectiilor.
1500 GOTO 1210!Reia ajustarea pasului.
1510 !
1520 ! ==========   P R O G R E S   I N  C A U T A R E   =========
1530 !
1540 ITN=ITN+1
1550 ALM=ALM/2
1560 INEF=0
1570 CALL CLEAR
1580 FOR I=1 TO NPAR
1590 W(I)=V(I)+AL*X(I)
1600 !LOCATE 13+I,1
1610 DISPLAY AT(13+I,1):USING "v(#)=###.### ":I,V(I)
1620 !LOCATE 13 + I, 65
1630 !PRINT USING " w()=."; I; w(I)
1640 NEXT I
1650 DISPLAY AT(1,1):USING "alfa=###.###^^^^":AL
1660 DISPLAY AT(2,1):USING "Par.L-M=##.##^^^^":ALM
1670 DISPLAY AT(3,1):USING "itn=#####":ITN
1680 DISPLAY AT(4,1):USING "rez=#####.##":F0
1690 NQ=0! Testeaza  cifrele semnif. in parametrii calculati.
1700 FOR I=1 TO NPAR
1710 IF ABS(V(I)/W(I)-1)<PREC THEN NQ=NQ+1
1720 NEXT I
1730 IF NQ=NPAR OR F0<SZ THEN 2060
1740 GOSUB 3270
1750 !
1760 !===========  CAUTAREA DE-A LUNGUL GRADIENTULUI  ============
1770 !
1780 !-------------------   O B S E R V A T I E ------------------
1790 !      In aceasta sectiune  urmarim  sa  obtinem  un  progres
1800 ! orict  de   mic  in   minimizarea   functionalei  c.m.m.p.
1810 !      Justificarea  procedeului  a  fost data de Kontorovici
1820 ! care a observat ca matricea Hessianului  nu se schimba prea
1830 ! mult in vecinatatea solutiei. Procedeul este in mod special
1840 ! util daca parametrul Levenberg-Marquart este nenul.
1850 !------------------------------------------------------------
1860 FOR I=1 TO NPAR
1870 !V(I)=W(I)
1880 W(I)=W(I)+X(I)*AL/2
1890 NEXT I
1900 GOSUB 2480
1910 IF F<F0 THEN F0=F:GOTO 1650
1920 !
1930 !==========  ACTUALIZAREA PARAMETRILOR OPTIMI ===============
1940 !
1950 FOR I=1 TO NPAR
1960 V(I)=VT(I)
1970 NEXT I
1980 KO=0
1990 INEF=0
2000 !===================  RELUAREA CAUTARII  ====================
2010 GOTO 780
2020 !
2030 !===================  SFIRSITUL CAUTARII  ===================
2040 !
2050 GOSUB 3200! Actualizarea param. optimi pt. raportul final.
2060 !FOR I=1 TO 7
2080 !NEXT I
2090 GOSUB 2930
2100 CALL CLEAR
2110 !LOCATE 1, 30
2120 DISPLAY AT(1,1):"TERMINAT"
2130 IF NQ<>NPAR THEN 2170
2140 A$="TERMINATIE NORMALA."
2150 A$=A$&" S-AU OBTINUT ## CIFRE SEMNIF. IN ### ITERATII"
2160 DISPLAY AT(2,1):USING A$:NSIG,ITN :: GOTO 2230
2170 A$="TERMINATIE ANORMALA "
2180 IF INEF=0 THEN 2210
2190 PRINT A$&"datorata cautarii ineficiente"
2200 GOTO 2230
2210 IF IFN<MAXFN THEN 2230
2220 PRINT A$&". S-a depasit nr. maxim de eval. ale fct."
2230 FOR I=1 TO NPAR
2240 Z$="V(#)=##.####^^^^+/-#.###^^^^"
2250 DISPLAY AT(5+I,1):USING Z$:I,W(I),SQR(ABS(AINV(I,I)))
2260 NEXT I
2270 !PRINT "Apasa orice tasta pentru vizualizarea finala"
2280 GOTO 2280
2290 CALL CLEAR
2300 GOSUB 3270
2310 PRINT "Apasa orice tasta"
2320 !CALL KEY() = "" THEN 2320
2330 END
2340 !************************************************************
2350 !*          SFIRSITUL  PROGRAMULUI  PRINCIPAL               *
2360 !************************************************************
2370 !
2380 !************************************************************
2390 !   SUBRUTINA PT. SCH. SEMN. IN CAZUL CAUTARII INEFICIENTE
2400 !************************************************************
2410 FOR I=1 TO NPAR
2420 X(I)=-X(I)/4
2430 NEXT I
2440 RETURN
2450 !************************************************************
2460 !      SUBRUTINA PENTRU CALCULUL FUNCTIONALEI C.M.M.P.
2470 !************************************************************
2480 F=0
2490 IFN=IFN+1
2500 GAMW=W(2)^2/4
2510 FOR I=1 TO NPCT
2520 AUX=W(1)-W(2)*W(3)/((I-W(4))^2+GAMW)
2530 AUX=AUX-W(3)*W(2)/((I-W(5)*2+W(4))^2+GAMW)
2540 DC(I)=AUX
2550 AUX=D(I)-AUX
2560 F=F+AUX*AUX
2570 NEXT I
2580 RETURN
2590 !************************************************************
2600 !       SUBRUTINA PENTRU REZOLVAREA SISTEMULUI LINIAR
2610 !************************************************************
2620 !--Procedeul ales este elim. Gaussiana cu pivotare partiala--
2630 !
2640 FOR I=1 TO NPAR! Initializare
2650 A(I,I)=A(I,I)+ALM
2660 !          Parametrul L.-M. este adaugat pe diag. principala.
2670 NEXT I
2680 FOR K=1 TO NPAR-1! Contorizarea eliminarilor.
2690 !
2700 !------------  Triunghiularizarea matricii A  --------------
2710 !
2720 FOR I=K+1 TO NPAR
2730 T=A(I,K)/A(K,K)
2740 FOR J=K+1 TO NPAR
2750 A(I,J)=A(I,J)-T*A(K,J)
2760 NEXT J
2770 B(I)=B(I)-T*B(K)
2780 NEXT I
2790 NEXT K
2800 !
2810 !--------  Construirea solutiilor sistemului liniar   -------
2820 !
2830 X(NPAR)=B(NPAR)/A(NPAR,NPAR)
2840 FOR I=NPAR-1 TO 1 STEP-1
2850 S=0
2860 FOR J=I+1 TO NPAR
2870 S=S+A(I,J)*X(J)
2880 NEXT J
2890 X(I)=(B(I)-S)/A(I,I)
2900 NEXT I
2910 RETURN
2920 !
2930 !************************************************************
2940 !     SUBRUTINA PENTRU CALCULUL INVERSEI MAT. HESSIAN
2950 !************************************************************
2960 FOR II=1 TO NPAR
2970 FOR J=1 TO NPAR
2980 B(J)=0
2990 NEXT J
3000 B(II)=1
3010 AINV(II,NPAR)=B(NPAR)/A(NPAR,NPAR)
3020 FOR I=NPAR-1 TO 1 STEP-1
3030 S=0
3040 FOR J=I+1 TO NPAR
3050 S=S+A(I,J)*AINV(II,J)
3060 NEXT J
3070 AINV(II,I)=(B(I)-S)/A(I,I)
3080 NEXT I
3090 NEXT II
3100 RETURN
3110 !************************************************************
3120 !      SUBRUTINA PT. INCERC. DIF. PASI IN OPTIM. PARAM.
3130 !************************************************************
3140 FOR I=1 TO NPAR
3150 W(I)=V(I)+ALFA*X(I)
3160 NEXT I
3170 GOSUB 2480
3180 RETURN
3190 !
3200 !************************************************************
3210 ! SUB. PT.CONSTRUIREA VECT. DE LUCRU CU VAL.FINALA A PASULUI
3220 !************************************************************
3230 FOR I=1 TO NPAR
3240 W(I)=V(I)+AL*X(I)
3250 NEXT I
3260 RETURN
3270 !************************************************************
3280 !  SUBRUTINA PENTRU REPREZENTAREA GRAFICA A SPECTRULUI
3290 !************************************************************
3300 ! FOR I = 1 TO NPCT
3310 !   PP2 = (D(I) - MI) / SCL * 150
3320 !   PSET (I * 5, 180 - PP2 + 1)
3330 !   PSET (I * 5, 180 - PP2 - 1)
3340 !   PSET (I * 5, 180 - PP2)
3350 !   PSET (I * 5 - 1, 180 - PP2)
3360 !   PSET (I * 5 + 1, 180 - PP2)
3370 ! NEXT I
3380 !PSET (1, 180 - (DC(1) - MI) / SCL * 150)
3390 !FOR I = 2 TO NPCT
3400 !  LINE -(I * 5, 1 * (180 - (DC(I) - MI) / SCL * 150))
3410 !NEXT I
3420 RETURN