src/snark/rtfort.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/rtfort.cpp $
   5  $LastChangedRevision: 95 $
   6  $Date: 2014-07-02 19:43:36 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  ************************************************************
  11
  12  SUBROUTINE RTFORT(X, N, INVDIR)
  13
  14  PURPOSE
  15  RTFORT COMPUTES THE HERMITIAN FOURIER SPECTRUM OF A ONE OR TWO
  16  DIMENSIONAL ARRAY OF REAL DATA (IF INVDIR=1), OR  THE REAL DATA
  17  CORRESPONDING TO THIS SPECTRUM (IF INVDIR=-1).
  18
  19  USAGE
  20  TO COMPUTE DIRECT FOURIER TRANSFORM OF REAL ARRAY X (INVDIR=1)
  21  CALL RTFORT( X, N, INVDIR )
  22  HALF OF THE CONJUGATE-SYMMETRIC (HERMITIAN) SPECTRUM
  23  IS NOW STORED IN X AS COMPLEX FOURIER COEFFICIENTS.
  24  TO COMPUTE INVERSE TRANSFORM OF SEMI-SPECTRUM (INVDIR = -1)
  25  CALL RTFORT(X, N, INVDIR)
  26  ARRAY X NOW CONTAINS REAL DATA.
  27
  28  PARAMETERS
  29  X
  30  ONE DIMENSIONAL, N(1)=1 : X CONTAINS 2*N(2) REAL NUMBERS IN
  31  CONSECUTIVE LOCATIONS AS INPUT TO DIRECT TRANSFORM OR AS RESULT
  32  OF INVERSE TRANSFORM.  THE SEMI-SPECTRUM CONSISTS OF ( N(2) + 1
  33  COMPLEX COEFFICIENTS STORED IN 2*N(2) + 2 LOCATIONS OF X.
  34
  35  TWO DIMENSIONAL : X IS A REAL ARRAY OF (FAST*SLOW) DIMENSIONS
  36  2N(2) * N(1) AS INPUT TO DIRECT TRANSFORM OR AS RESULT OF
  37  INVERSE TRANSFORM.  THE SEMI-SPECTRUM CONSISTS OF
  38  ( N(2) + 1 ) * N(1)  (FAST*SLOW) COMPLEX FOURIER COEFFICIENTS
  39  STORED IN ( 2*N(2) + 2 ) * N(1) LOCATIONS OF X.
  40
  41  N.B.  X MUST BE DIMENSIONED TO ACCOMODATE THE SEMI-SPECTRUM,
  42  WHICH REQUIRES MORE STORAGE THAN ITS CORRESPONDING REAL ARRAY.
  43
  44  N
  45  A 3-ELEMENT VECTOR, N(1)= SLOW DIMENSION OF X ARRAY
  46  2*N(2) = FAST DIMENSION OF REAL X ARRAY (AN EVEN NUMBER)
  47  N(3) = 1 (ONLY 2-DIMENSIONAL ARRAYS ALLOWED)
  48  THE VALUES OF N(1) AND N(2) MUST BE COMPATIBLE WITH THE
  49  PARTICULAR COMPLEX FFT ROUTINE USED (E.G. SNFFT, HARM).
  50
  51  INVDIR = 1   FOR DIRECT (FORWARD) TRANSFORM
  52  = -1  FOR INVERSE TRANSFORM
  53
  54  SUBROUTINE CALLED:  TRANSM
  55
  56  BASED ON R2FORT SUBROUTINE WRITTEN BY T.M.PETERS, JUNE 1971,
  57  FOR USE WITH RADIX TWO FFT.  RTFORT EXTENDS CAPABILITY TO
  58  RADIX ODD TRANSFORMS.        R.M. LEWITT, SEPT 1977.
  59
  60  **************************************************************
  61  */
  62
  63 #include <cstdio>
  64 #include <cmath>
  65
  66 #include "blkdta.h"
  67 #include "snfft.h"
  68 #include "transm.h"
  69 #include "consts.h"
  70
  71 #include "rtfort.h"
  72
  73 void rtfort(REAL* x, INTEGER* nin, INTEGER invdir)
  74 {
  75
  76         INTEGER n[3];
  77
  78         INTEGER m;
  79         INTEGER ns;
  80         INTEGER m1;
  81         INTEGER mm;
  82         INTEGER mm2;
  83         REAL sd;
  84         REAL cd;
  85         REAL sn;
  86         REAL cn;
  87         INTEGER if_1;
  88         INTEGER ib;
  89         REAL ct;
  90         REAL ca;
  91         INTEGER mtype;
  92         INTEGER index;
  93         INTEGER k;
  94         INTEGER ishift;
  95         INTEGER i;
  96         INTEGER igap;
  97         INTEGER jt;
  98         INTEGER i1;
  99         INTEGER j1;
 100         INTEGER j2;
 101         INTEGER k1;
 102         INTEGER nl;
 103
 104         n[0] = nin[0];
 105         n[1] = nin[1];
 106         n[2] = 1;
 107
 108         if (invdir == 1)
 109         {
 110                 snfft(x, n, invdir);
 111         }
 112
 113         m = n[1];
 114         ns = n[0];
 115         m1 = m + 1;
 116         mm = m + m;
 117         mm2 = mm + 2;
 118
 119         // INITIALISE PARAMETERS FOR RECURSIVE GENERATION OF SINES
 120         // AND COSINES
 121         sd = Consts.pid2 / m;
 122         cd = ((REAL) 2.0) * SQR((REAL) sin(sd));
 123
 124         sd = (REAL) -sin(sd + sd);
 125
 126         // ASSUMING THE DIRECT FOURIER TRANSFORM HAS EXP(-I...)
 127         // IF DIRECT FT ROUTINE HAS EXP(+I...) FACTOR, THEN
 128         // REPLACE ABOVE STATEMENT BY
 129         // SD=SIN(SD+SD)
 130
 131         sn = 0.;
 132         if (!(ns > 1))
 133         {
 134
 135                 // ONE DIMENSIONAL CASE
 136                 if (!(invdir < 0))
 137                 {
 138
 139                         // FORWARD TRANSFORM
 140                         cn = 1.;
 141
 142                         // REPEAT FIRST TWO REAL NUMBERS AT ARRAY END FOR PERIODICITY.
 143                         x[mm + 0] = x[0];
 144                         x[mm + 1] = x[1];
 145                 }
 146                 else
 147                 {
 148
 149                         // INVERSE TRANSFORM
 150                         cn = -1.;
 151                         sd = -sd;
 152                 }
 153
 154                 // APPLY TRANSFORM SYMMETRY RELATIONS
 155                 for (if_1 = 0; if_1 < m1; if_1 += 2)
 156                 {
 157                         ib = mm - if_1;
 158                         transm(x, if_1, ib, sn, cn);
 159                         ct = cn - (cd * cn + sd * sn);
 160                         sn = sn + (sd * cn - cd * sn);
 161                         ca = ((REAL) 0.5) / (SQR(ct) + SQR(sn)) + ((REAL) 0.5);
 162                         sn = ca * sn;
 163                         cn = ca * ct;
 164                 }
 165
 166                 if (!(invdir > 0))
 167                 {
 168
 169                         // FOR INVERSE TRANSFORM, ZERO LAST TWO ARRAY ELEMENTS
 170                         x[mm + 0] = 0.;
 171                         x[mm + 1] = 0.;
 172                 }
 173         }
 174         else
 175         {
 176                 // TWO DIMENSIONAL CASE
 177                 // MTYPE = 0 (1)  FOR  M ODD (EVEN)
 178                 mtype = (m / 2) * 2 - m + 1;
 179                 if (!(invdir < 0))
 180                 {
 181                         cn = 1.;
 182
 183                         // FORWARD TRANSFORM. EXPAND ARRAY TO TRANSFORM FORMAT
 184                         index = ns * mm;
 185                         for (k = 1; k <= ns; k++)
 186                         {
 187                                 ishift = (ns - k) * 2;
 188                                 for (i = 0; i < mm; i++)
 189                                 {
 190                                         index--;
 191                                         x[index + ishift] = x[index];
 192                                 }
 193                                 igap = index + ishift + mm;
 194                                 x[igap] = x[index];
 195                                 x[igap + 1] = x[index + 1];
 196                         }
 197                 }
 198                 else
 199                 {
 200
 201                         // INVERSE TRANSFORM
 202                         cn = -1.;
 203                         sd = -sd;
 204                 }
 205
 206                 jt = (ns + 1) * mm2;
 207                 for (i1 = 0; i1 < m; i1 += 2)
 208                 {
 209                         // K=0 CASE SEPARATELY
 210                         j1 = i1;
 211                         j2 = mm - i1;
 212                         transm(x, j1, j2, sn, cn);
 213
 214                         index = 0;
 215                         for (k1 = 1; k1 < ns; k1++)
 216                         {
 217                                 index += mm2;
 218                                 j1 = index + i1;
 219                                 j2 = jt - j1 - 2;
 220                                 transm(x, j1, j2, sn, cn);
 221                         }
 222
 223                         ct = cn - (cd * cn + sd * sn);
 224                         sn = sn + (sd * cn - cd * sn);
 225                         ca = ((REAL) 0.5) / (SQR(ct) + SQR(sn)) + ((REAL) 0.5);
 226                         sn = ca * sn;
 227                         cn = ca * ct;
 228                 }
 229
 230                 if (!(mtype == 0))
 231                 {
 232
 233                         // M EVEN, TAKE I1=M1 AS SPECIAL CASE
 234                         m1 -= 1;             // decrement m1 - swr 11/25/04
 235                         i1 = m1;
 236                         transm(x, m1, m1, sn, cn);
 237                         nl = ns / 2 + 1;
 238                         index = 0;
 239                         for (k1 = 2; k1 <= nl; k1++)
 240                         {   // k => k1 - swr 11/25/04
 241                                 index += mm2;
 242                                 j1 = index + m1;
 243                                 j2 = jt - j1 - 2;  // decrement by 2 - swr 11/25/04
 244                                 transm(x, j1, j2, sn, cn);
 245                         }
 246                 }
 247
 248                 if (!(invdir > 0))
 249                 {
 250
 251                         // INVERSE TRANSFORM, REVERT TO ORIGINAL FORMAT
 252                         index = mm;
 253                         for (k = 1; k < ns; k++)
 254                         {
 255                                 ishift = k * 2;
 256                                 for (i = 1; i <= mm; i++)
 257                                 {
 258                                         //index++;
 259                                         x[index] = x[index + ishift];
 260                                         index++;
 261                                 }
 262                         }
 263
 264                         // ZERO REMAINDER OF REAL ARRAY
 265                         j1 = index;
 266                         j2 = ns * mm2;
 267                         for (i = j1; i < j2; i++)
 268                         {
 269                                 x[i] = 0.;
 270                         }
 271                 }
 272         }
 273
 274         if (invdir == -1)
 275         {
 276                 snfft(x, n, invdir);
 277         }
 278
 279         return;
 280 }