From 5ca9582e06e0d697003df6a41d0c479fc2bbd461 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Fri, 15 Dec 2023 14:20:26 +0100
Subject: [PATCH] Implement status logging in FORTRAN code

---
 src/cluster/clu.f     | 31 ++-------------------
 src/sphere/sph.f      | 63 ++-----------------------------------------
 src/sphere/sphere.cpp |  2 +-
 3 files changed, 5 insertions(+), 91 deletions(-)

diff --git a/src/cluster/clu.f b/src/cluster/clu.f
index 9c902231..2fb83995 100644
--- a/src/cluster/clu.f
+++ b/src/cluster/clu.f
@@ -330,12 +330,8 @@ C     1", arg =",ARG
       GO TO 490
  132  CONTINUE
       CALL CMS(AM)
-      CALL SUMMAT(AM,CCSAM)
-      PRINT *,"DEBUG: after CMS CCSAM =",CCSAM
       NDIT=2*NSPH*NLIM
       CALL LUCIN(AM,MXNDM,NDIT,JER)
-      CALL SUMMAT(AM,CCSAM)
-      PRINT *,"DEBUG: after LUCIN CCSAM =",CCSAM
       IF(JER.EQ.1)GO TO 495
       CALL ZTM(AM)
 
@@ -358,7 +354,6 @@ C
   158 WRITE(IW,6061)
   160 CS0=0.25D0*VK*VK*VK/DACOS(C0)
       CALL SCR0(VK,EXRI)
-      PRINT *,"DEBUG: after SCR0 TFSAS =",TFSAS
       SQK=VK*VK*EXDC
       CALL APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS)
       CALL RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -639,6 +634,7 @@ C
   482 IF(ISAM.LE.1)THSL=THSL+THSSTP
   484 PH=PH+PHSTP
   486 TH=TH+THSTP
+      PRINT *,"INFO: done scale iteration ",JXI
   488 CONTINUE
       GO TO 500
   490 WRITE(IW,6550)JER,LCALC,ARG
@@ -650,6 +646,7 @@ C
       CLOSE(IW)
       CLOSE(IT)
       CLOSE(ITIN)
+      PRINT *,"INFO: output written to OCLU."
       STOP
       END
       SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -3504,28 +3501,4 @@ CCC   DIMENSION AV(NDDMST*NDDMST)
  3    CDTP=CDTP+AV((J-1)*ISTEP+I)*AV(J+IR)
       RETURN
       END
-      SUBROUTINE SUMMAT(AM,RSLT)
-      DIMENSION AM(921600)
-      COMPLEX*16 AM,CC0,RSLT
-      CC0=(0.0D0,0.0D0)
-      RSLT=CC0
-      DO 4 I=1,921600
- 4    RSLT=RSLT+AM(I)
-      RETURN
-      END
-      SUBROUTINE LOGMAT(AM)
-      DIMENSION AM(960,960)
-      COMPLEX*16 AM
-      IL=36
-      OPEN(IL,FILE="matrix.log",STATUS="UNKNOWN")
-      DO 5 I=1,960
-      DO 6 J=1,960
-      WRITE(IL,9901)I,J,DREAL(AM(I,J))
- 6    WRITE(IL,9902)I,J,DIMAG(AM(I,J))
- 5    CONTINUE
-      CLOSE(IL)
- 9901 FORMAT("R:",1I3,",",1I3,",",1D15.7)
- 9902 FORMAT("I:",1I3,",",1I3,",",1D15.7)
-      RETURN
-      END
 CCC
diff --git a/src/sphere/sph.f b/src/sphere/sph.f
index bebacd74..40bd15a5 100644
--- a/src/sphere/sph.f
+++ b/src/sphere/sph.f
@@ -270,7 +270,6 @@ C
 
       CS0=0.25D0*VK*VK*VK/PIGH
       CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI)
-      PRINT *,"DEBUG: TFSAS =", TFSAS
       SQK=VK*VK*EXDC
       CALL APS(ZPV,LM,NSPH,IOG,RMI,REI,SQK,GAPS)
       CALL RABAS(INPOL,LM,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -397,6 +396,7 @@ C
   482 IF(ISAM.LE.1)THSL=THSL+THSSTP
   484 PH=PH+PHSTP
   486 TH=TH+THSTP
+      PRINT *,"INFO: done scale iteration ",JXI
   488 CONTINUE
       GO TO 500
   490 WRITE(IW,6550)JER,LCALC,ARG
@@ -406,6 +406,7 @@ C
       CLOSE(IW)
       CLOSE(IT)
       CLOSE(ITIN)
+      PRINT *,"INFO: output written to OSPH."
       STOP
       END
       SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -456,19 +457,11 @@ CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
       TQSPE(2,I)=TQSPE(2,I)*PIG2
       TQSPS(1,I)=TQSPS(1,I)*PIG2
       TQSPS(2,I)=TQSPS(2,I)*PIG2
-      PRINT *,"DEBUG: TQSPE(1,",I,") =",TQSPE(1,I)
-      PRINT *,"DEBUG: TQSPE(2,",I,") =",TQSPE(2,I)
-      PRINT *,"DEBUG: TQSPS(1,",I,") =",TQSPS(1,I)
-      PRINT *,"DEBUG: TQSPS(2,",I,") =",TQSPS(2,I)
       GO TO 80
    75 TQSE(1,I)=TQSE(1,I)*PIG2
       TQSE(2,I)=TQSE(2,I)*PIG2
       TQSS(1,I)=TQSS(1,I)*PIG2
       TQSS(2,I)=TQSS(2,I)*PIG2
-      PRINT *,"DEBUG: TQSE(1,",I,") =",TQSE(1,I)
-      PRINT *,"DEBUG: TQSE(2,",I,") =",TQSE(2,I)
-      PRINT *,"DEBUG: TQSS(1,",I,") =",TQSS(1,I)
-      PRINT *,"DEBUG: TQSS(2,",I,") =",TQSS(2,I)
    80 CONTINUE
       RETURN
       END
@@ -584,7 +577,6 @@ CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
       PIGFSQ=6.4D+1*DACOS(C0)**2
       CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
       NLMM=LM*(LM+2)
-C     PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1)
       DO 14 I=1,NSPH
       IOGI=IOG(I)
       IF(IOGI.LT.I)GO TO 14
@@ -596,30 +588,14 @@ C     PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1)
       DO 10 L=1,LM
       RM=1.0D0/RMI(L,I)
       RE=1.0D0/REI(L,I)
-C     PRINT *,"DEBUG: RM =",RM
-C     PRINT *,"DEBUG: RE =",RE
       LTPO=L+L+1
       DO 10 IM=1,LTPO
       K=K+1
       KE=K+NLMM
-C     PRINT *,"DEBUG: W(",K,",3) =",W(K,3)
-C     PRINT *,"DEBUG: W(",K,",1) =",W(K,1)
-C     PRINT *,"DEBUG: W(",KE,",3) =",W(KE,3)
-C     PRINT *,"DEBUG: W(",KE,",1) =",W(KE,1)
       S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE
-C     PRINT *,"DEBUG: S11 =",S11
-C     PRINT *,"DEBUG: W(",K,",4) =",W(K,4)
-C     PRINT *,"DEBUG: W(",KE,",4) =",W(KE,4)
       S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE
-C     PRINT *,"DEBUG: W(",K,",2) =",W(K,2)
-C     PRINT *,"DEBUG: W(",KE,",2) =",W(KE,2)
       S12=S12-W(K,3)*W(K,2)*RM-W(KE,3)*W(KE,2)*RE
    10 S22=S22-W(K,4)*W(K,2)*RM-W(KE,4)*W(KE,2)*RE
-C     PRINT *,"DEBUG: CSAM =",CSAM
-C     PRINT *,"DEBUG: S11 =",S11
-C     PRINT *,"DEBUG: S21 =",S21
-C     PRINT *,"DEBUG: S12 =",S12
-C     PRINT *,"DEBUG: S22 =",S22
       SAS(I,1,1)=S11*CSAM
       SAS(I,2,1)=S21*CSAM
       SAS(I,1,2)=S12*CSAM
@@ -670,7 +646,6 @@ CCC  1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH)
      1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL
    30 CONTINUE
       GAPS(I)=DREAL(SUM)*COFS
-      PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I)
    40 CONTINUE
       RETURN
       END
@@ -781,10 +756,6 @@ CCC   CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M)
       SINT=DSIN(TH)
       COSP=DCOS(PH)
       SINP=DSIN(PH)
-C     PRINT *,"DEBUG: cost =",COST
-C     PRINT *,"DEBUG: sint =",SINT
-C     PRINT *,"DEBUG: cosp =",COSP
-C     PRINT *,"DEBUG: sinp =",SINP
       U(1)=COSP*SINT
       U(2)=SINP*SINT
       U(3)=COST
@@ -828,8 +799,6 @@ CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
       DO 60 N=1,NSPH
    60 ARG(N)=-ARG(N)
  65   CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
-C     PRINT *,"DEBUG: in WMAMP and calling PWMA with lm =",LM,
-C    1"and iis =",IIS
       CALL PWMA(UP,UN,YLM,INPOL,LM,IIS)
       RETURN
       END
@@ -928,11 +897,9 @@ CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
    55 ARGS(N)=-COSTS*RZZ(N)
    60 CONTINUE
  75   CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
-C     PRINT *,"DEBUG: in WMASP and calling PWMA"
       CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ)
       IF(IBF.LT.0)RETURN
       CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM)
-C     PRINT *,"DEBUG: in WMASP and calling PWMA with IBF =",IBF
       CALL PWMA(UPS,UNS,YLM,INPOL,LM,2)
       RETURN
       END
@@ -963,7 +930,6 @@ CCC   DIMENSION YLM(NLWM+2)
       CM2=.5D0*DCMPLX(UN(1),UN(2))
       CP2=.5D0*DCMPLX(UN(1),-UN(2))
       CZ2=UN(3)
-C     PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2)
       DO 20 L=1,LW
       LF=L+1
       LFTL=LF*L
@@ -980,35 +946,26 @@ C     PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2)
       CM=DSQRT(X)
       CZ=M
       W(K,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
-C     IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
  20   W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL
-C 20  PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
       DO 30 K=1,NLWM
       I=K+NLWM
       W(I,ISPO)=UIM*W(K,ISPT)
-C     PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
  30   W(I,ISPT)=-UIM*W(K,ISPO)
-C 30  PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
       IF(INPOL.EQ.0)GO TO 42
       DO 40 K=1,NLWM
       I=K+NLWM
       C1=(W(K,ISPO)+UIM*W(K,ISPT))*SQRTWI
       C2=(W(K,ISPO)-UIM*W(K,ISPT))*SQRTWI
       W(K,ISPO)=C2
-C     PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
       W(I,ISPO)=-C2
-C     PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
       W(K,ISPT)=C1
-C     PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
  40   W(I,ISPT)=C1
-C 40  PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
    42 IF(ISQ.EQ.0)RETURN
       DO 50 I=1,2
       IPT=I+2
       IPIS=I+IS
       DO 50 K=1,NLWMT
  50   W(K,IPT)=DCONJG(W(K,IPIS))
-C 50  PRINT *,"DEBUG: W(",K,",",IPT,") =",W(K,IPT)
       RETURN
       END
       SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH)
@@ -1101,9 +1058,7 @@ CCC  2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
       CCNC=FBI(LPO)*DFB
       CCND=FB(LPO)*DFBI
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
-C     PRINT *,"DEBUG: gone 60, RMI(",L,",",I,") =",RMI(L,I)
  60   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
-C 60  PRINT *,"DEBUG: gone 60, REI(",L,",",I,") =",REI(L,I)
       RETURN
    65 DO 80 L=1,LI
       LPO=L+1
@@ -1138,7 +1093,6 @@ cccccccccccccccccccccc
    80 DREF(L)=DY2*SZ+Y2
       CRI=(1.0D0,0.0D0)
       IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC
-C     PRINT *,"DEBUG: going 90, AREX =",AREX
       DO 90 L=1,LI
       LPO=L+1
       LTPO=LPO+L
@@ -1150,13 +1104,11 @@ C     PRINT *,"DEBUG: going 90, AREX =",AREX
       CCNC=RMF(L)*DFB
       CCND=DRMF(L)*FB(LPO)*SZ*LTPO
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
-C      PRINT *,"DEBUG: gone 90, RMI(",L,",",I,") =",RMI(L,I)
       CCNA=REF(L)*DFN
       CCNB=DREF(L)*FN(LPO)*SZ*LTPO
       CCNC=REF(L)*DFB
       CCND=DREF(L)*FB(LPO)*SZ*LTPO
  90   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
-C 90  PRINT *,"DEBUG: gone 90, REI(",L,",",I,") =",REI(L,I)
       RETURN
       END
       SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2)
@@ -1498,18 +1450,13 @@ CCC   LL=LM
       PI4=DACOS(0.0D0)*8.0D0
       PI4IRS=1.0D0/DSQRT(PI4)
       X=COSRTH
-C     PRINT *,"DEBUG: X =",X
       Y=DABS(SINRTH)
-C     PRINT *,"DEBUG: Y =",Y
       CLLMO=3.0D0
       CLL=1.5D0
       YTOL=Y
       PLEGN(1)=1.0D0
-C     PRINT *,"DEBUG: PLEGN( 1 ) =",PLEGN(1)
       PLEGN(2)=X*DSQRT(CLLMO)
-C     PRINT *,"DEBUG: PLEGN( 2 ) =",PLEGN(2)
       PLEGN(3)=YTOL*DSQRT(CLL)
-C     PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3)
       SINRMP(1)=SINRPH
       COSRMP(1)=COSRPH
       IF(LL.LT.2)GO TO 30
@@ -1529,17 +1476,14 @@ C     PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3)
       CDM=(LTMO-2)*LS
  10   PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
      1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)
-C10   PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK)
       LPK=L+K
       CLTPO=LTPO
       PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO)
-C     PRINT *,"DEBUG: PLEGN(",LPK,") =",PLEGN(LPK)
       K=LPK+1
       CLT=LTPO-1
       CLL=CLL*(CLTPO/CLT)
       YTOL=YTOL*Y
       PLEGN(K)=YTOL*DSQRT(CLL)
-C     PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K)
       SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO)
    20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO)
    30 L=0
@@ -1548,17 +1492,14 @@ C     PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K)
       L0Y=K+1
       L0P=K/2+1
       YLM(L0Y)=PI4IRS*PLEGN(L0P)
-C     PRINT *, "DEBUG: YLM(",L0Y,") =",YLM(L0Y)
       GO TO 45
    44 LMP=L0P+M
       SAVE=PI4IRS*PLEGN(LMP)
       LMY=L0Y+M
       YLM(LMY)=SAVE*DCMPLX(COSRMP(M),SINRMP(M))
       IF(MOD(M,2).NE.0)YLM(LMY)=-YLM(LMY)
-C     PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
       LMY=L0Y-M
       YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M))
-C     PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
  45   IF(M.GE.L)GO TO 47
       M=M+1
       GO TO 44
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index e5adf042..2cbafe99 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -20,7 +20,7 @@ using namespace std;
 void sphere(string config_file, string data_file, string output_path) {
 	complex<double> arg, s0, tfsas;
 	double th, ph;
-	printf("INFO: making legacy configuration ...\n");
+	printf("INFO: making legacy configuration...\n");
 	ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
 	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
 	if (sconf->number_of_spheres == gconf->number_of_spheres) {
-- 
GitLab