SWIN <- function(nn,fsin,df,band){
    nfold <- nn/2+1
    tt    <- 1.0/df
    udf   <- 280.0/151.0/band*df
    lmax  <- floor(2.0/udf)+1
    dif   <- numeric(lmax-1)
    w     <- numeric(lmax)
    g     <- numeric(nfold)
    fs    <- numeric(nfold)
    fs    <- fsin[1:nfold]

    if(udf<=0.5){
        ii            <- seq(1,lmax-1,by=1)
        dif[1:lmax-1] <- 0.5*pi*ii*udf
        w[1]          <- 0.75*udf
        w[2:lmax]     <- w[1]*(sin(dif)/dif)^4.0
        g[1:nfold]    <- 2.0*fs[1:nfold]^2/tt
        g[1]          <- fs[1]^2/tt
        g[nfold]      <- fs[nfold]^2/tt

        ll <- lmax*2-1
        ln <- ll-1+nfold
        lt <- (ll-1)*2+nfold
        le <- lt-lmax+1
        g1 <- numeric(lt)
        g2 <- numeric(le)
        for(i in 1:nfold){
            wg         <- g[i]
            g1[ll-1+i] <- wg
        }
        for(i in lmax:le){
            s <- w[1]*g1[i]
            for(j in 2:lmax){
                s <- s+w[j]*(g1[i-j+1]+g1[i+j-1])
            }
            g2[i] <- s
        }
        for(j in 2:lmax){
            g2[ll+j-1] <- g2[ll+j-1]+g2[ll-j+1]
            g2[ln-j+1] <- g2[ln-j+1]+g2[ln+j-1]
        }
        for(i in 1:nfold){
            wg   <-g2[ll-1+i]
            g[i] <- wg
        }
        fs[1:nfold]        <- sqrt(g[1:nfold]*tt/2.0)
        fs[1]              <- sqrt(g[1]*tt)
        fs[nfold]          <- sqrt(g[nfold]*tt)
    }
    return(fs)
}
#******************
# main
#******************
for(iii in 1:3){
    if(iii==1)fnameR="dat_acc_TCGH16_EW2.csv"
    if(iii==2)fnameR="dat_acc_TCGH16_NS2.csv"
    if(iii==3)fnameR="dat_acc_TCGH16_UD2.csv"
    if(iii==1)fnameW="dat_fsp_1.txt"
    if(iii==2)fnameW="dat_fsp_2.txt"
    if(iii==3)fnameW="dat_fsp_3.txt"
    if(iii==1)strn="TCGH16_EW2"
    if(iii==2)strn="TCGH16_NS2"
    if(iii==3)strn="TCGH16_UD2"

    d  <- scan(fnameR,skip=3)
    k  <- length(d)
    nn <- 2
    while(nn<k)nn <- nn*2
    x      <- numeric(nn)
    x[1:k] <- d[1:k]
    z      <- fft(x)
    z1     <- Re(z)/nn
    z2     <- Im(z)/nn
    nfold  <- nn/2+1
    dt     <- 0.01
    band   <- 1.0
    df     <- 1.0/dt/nn
    fs1    <- numeric(nfold)
    fs2    <- numeric(nfold)
    frq    <- numeric(nfold)
    fs1    <- dt*nn*sqrt(z1[1:nfold]^2+z2[1:nfold]^2)
    fs2    <- SWIN(nn,fs1,df,band)
    ii     <- seq(0,nfold-1,by=1)
    frq    <- ii*df
    fsp    <- list(frq,fs1,fs2)
    names(fsp) <- c("frq","fs1","fsp")
    write.table(fsp,fnameW,quote=F,append=F,row.names=F)
}
#******************
# Plot
#******************
plot(
    0,0.1,
    log="y",
    type="n",
    xlim=c(0,20),
    ylim=c(0.1,10000),
    xlab="Frequency (Hz)",
    ylab="Fourier spectrum (gal * sec)",
    axes=FALSE
)
cc<-seq(1,9,by=1)
c1<-c(0,5,10,15,20)
c2<-c(cc*0.1,cc*1,cc*10,cc*100,cc*1000,10000)
abline(v=c1,col="gray")
abline(h=c2,col="gray")
axis(1,c(0,5,10,15,20)) 
axis(2,c(0.1,1,10,100,1000,10000),c("0.1","1","10","100","1000","10000")) 
box()
lines(frq,fs1,lwd=1,col="red")
lines(frq,fs2,lwd=2,col="black")
text(18,10000,strn)