/* Eddy Reconstruction with Exponential Outer Decay Fitting */
/* For Iron Spring PL/I Compiler */
 
EDDY_RECONSTRUCT: PROCEDURE OPTIONS(MAIN);
 
   DCL MAXN     FIXED BIN(31) STATIC INIT(1000);
   DCL X(MAXN), Y(MAXN), VT(MAXN) FLOAT DECIMAL(15);
   DCL N        FIXED BIN(31);
   DCL I        FIXED BIN(31);
   DCL XC, YC   FLOAT DECIMAL(15);
   DCL VT_FIT(MAXN) FLOAT DECIMAL(15);
   DCL LAMBDA   FLOAT DECIMAL(15);
 
   CALL GEN_SAMPLE_DATA(X, Y, VT, N);
 
   XC = 0.0; YC = 0.0;
   CALL FIT_OUTER_EXP(X, Y, VT, N, XC, YC, LAMBDA, VT_FIT);
 
   PUT SKIP LIST('Exponential decay fit λ =', LAMBDA);
   DO I = 1 TO N;
      PUT SKIP EDIT(X(I), Y(I), VT(I), VT_FIT(I)) 
         (F(10,4), X, F(10,4), X, F(10,4), X, F(10,4));
   END;
 
END EDDY_RECONSTRUCT;
 
FIT_OUTER_EXP: PROCEDURE(X, Y, VT, N, XC, YC, LAMBDA, VT_FIT);
   DCL X(*), Y(*), VT(*), VT_FIT(*) FLOAT DECIMAL(15);
   DCL N FIXED BIN(31);
   DCL XC, YC, LAMBDA FLOAT DECIMAL(15);
 
   DCL I FIXED BIN(31);
   DCL DX, DY, R, SUM_R, SUM_LOG FLOAT DECIMAL(15);
 
   SUM_R = 0.0; SUM_LOG = 0.0;
   DO I = 1 TO N;
      DX = X(I) - XC;
      DY = Y(I) - YC;
      R = SQRT(DX*DX + DY*DY);
      IF VT(I) > 0.0 THEN DO;
         SUM_R = SUM_R + R;
         SUM_LOG = SUM_LOG + LOG(VT(I));
      END;
   END;
 
   IF SUM_R > 0.0 THEN
      LAMBDA = -SUM_LOG / SUM_R;
   ELSE
      LAMBDA = 0.0;
 
   DO I = 1 TO N;
      DX = X(I) - XC;
      DY = Y(I) - YC;
      R = SQRT(DX*DX + DY*DY);
      VT_FIT(I) = EXP(-LAMBDA * R);
   END;
 
END FIT_OUTER_EXP;
 
GEN_SAMPLE_DATA: PROCEDURE(X, Y, VT, N);
   DCL X(*), Y(*), VT(*) FLOAT DECIMAL(15);
   DCL N FIXED BIN(31);
 
   DCL I FIXED BIN(31);
 
   N = 10;
   DO I = 1 TO N;
      X(I) = FLOAT(I - 5);
      Y(I) = 0.0;
      VT(I) = EXP(-ABS(X(I)));
   END;
 
END GEN_SAMPLE_DATA;
				LyogRWRkeSBSZWNvbnN0cnVjdGlvbiB3aXRoIEV4cG9uZW50aWFsIE91dGVyIERlY2F5IEZpdHRpbmcgKi8KLyogRm9yIElyb24gU3ByaW5nIFBML0kgQ29tcGlsZXIgKi8KCkVERFlfUkVDT05TVFJVQ1Q6IFBST0NFRFVSRSBPUFRJT05TKE1BSU4pOwoKICAgRENMIE1BWE4gICAgIEZJWEVEIEJJTigzMSkgU1RBVElDIElOSVQoMTAwMCk7CiAgIERDTCBYKE1BWE4pLCBZKE1BWE4pLCBWVChNQVhOKSBGTE9BVCBERUNJTUFMKDE1KTsKICAgRENMIE4gICAgICAgIEZJWEVEIEJJTigzMSk7CiAgIERDTCBJICAgICAgICBGSVhFRCBCSU4oMzEpOwogICBEQ0wgWEMsIFlDICAgRkxPQVQgREVDSU1BTCgxNSk7CiAgIERDTCBWVF9GSVQoTUFYTikgRkxPQVQgREVDSU1BTCgxNSk7CiAgIERDTCBMQU1CREEgICBGTE9BVCBERUNJTUFMKDE1KTsKCiAgIENBTEwgR0VOX1NBTVBMRV9EQVRBKFgsIFksIFZULCBOKTsKCiAgIFhDID0gMC4wOyBZQyA9IDAuMDsKICAgQ0FMTCBGSVRfT1VURVJfRVhQKFgsIFksIFZULCBOLCBYQywgWUMsIExBTUJEQSwgVlRfRklUKTsKCiAgIFBVVCBTS0lQIExJU1QoJ0V4cG9uZW50aWFsIGRlY2F5IGZpdCDOuyA9JywgTEFNQkRBKTsKICAgRE8gSSA9IDEgVE8gTjsKICAgICAgUFVUIFNLSVAgRURJVChYKEkpLCBZKEkpLCBWVChJKSwgVlRfRklUKEkpKSAKICAgICAgICAgKEYoMTAsNCksIFgsIEYoMTAsNCksIFgsIEYoMTAsNCksIFgsIEYoMTAsNCkpOwogICBFTkQ7CgpFTkQgRUREWV9SRUNPTlNUUlVDVDsKCkZJVF9PVVRFUl9FWFA6IFBST0NFRFVSRShYLCBZLCBWVCwgTiwgWEMsIFlDLCBMQU1CREEsIFZUX0ZJVCk7CiAgIERDTCBYKCopLCBZKCopLCBWVCgqKSwgVlRfRklUKCopIEZMT0FUIERFQ0lNQUwoMTUpOwogICBEQ0wgTiBGSVhFRCBCSU4oMzEpOwogICBEQ0wgWEMsIFlDLCBMQU1CREEgRkxPQVQgREVDSU1BTCgxNSk7CgogICBEQ0wgSSBGSVhFRCBCSU4oMzEpOwogICBEQ0wgRFgsIERZLCBSLCBTVU1fUiwgU1VNX0xPRyBGTE9BVCBERUNJTUFMKDE1KTsKCiAgIFNVTV9SID0gMC4wOyBTVU1fTE9HID0gMC4wOwogICBETyBJID0gMSBUTyBOOwogICAgICBEWCA9IFgoSSkgLSBYQzsKICAgICAgRFkgPSBZKEkpIC0gWUM7CiAgICAgIFIgPSBTUVJUKERYKkRYICsgRFkqRFkpOwogICAgICBJRiBWVChJKSA+IDAuMCBUSEVOIERPOwogICAgICAgICBTVU1fUiA9IFNVTV9SICsgUjsKICAgICAgICAgU1VNX0xPRyA9IFNVTV9MT0cgKyBMT0coVlQoSSkpOwogICAgICBFTkQ7CiAgIEVORDsKCiAgIElGIFNVTV9SID4gMC4wIFRIRU4KICAgICAgTEFNQkRBID0gLVNVTV9MT0cgLyBTVU1fUjsKICAgRUxTRQogICAgICBMQU1CREEgPSAwLjA7CgogICBETyBJID0gMSBUTyBOOwogICAgICBEWCA9IFgoSSkgLSBYQzsKICAgICAgRFkgPSBZKEkpIC0gWUM7CiAgICAgIFIgPSBTUVJUKERYKkRYICsgRFkqRFkpOwogICAgICBWVF9GSVQoSSkgPSBFWFAoLUxBTUJEQSAqIFIpOwogICBFTkQ7CgpFTkQgRklUX09VVEVSX0VYUDsKCkdFTl9TQU1QTEVfREFUQTogUFJPQ0VEVVJFKFgsIFksIFZULCBOKTsKICAgRENMIFgoKiksIFkoKiksIFZUKCopIEZMT0FUIERFQ0lNQUwoMTUpOwogICBEQ0wgTiBGSVhFRCBCSU4oMzEpOwoKICAgRENMIEkgRklYRUQgQklOKDMxKTsKCiAgIE4gPSAxMDsKICAgRE8gSSA9IDEgVE8gTjsKICAgICAgWChJKSA9IEZMT0FUKEkgLSA1KTsKICAgICAgWShJKSA9IDAuMDsKICAgICAgVlQoSSkgPSBFWFAoLUFCUyhYKEkpKSk7CiAgIEVORDsKCkVORCBHRU5fU0FNUExFX0RBVEE7