wavelet.c 8.7 KB
Newer Older
tilman.metz's avatar
tilman.metz committed
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
tilman.metz's avatar
tilman.metz committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
tilman.metz's avatar
tilman.metz committed
5
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
tilman.metz's avatar
tilman.metz committed
7 8 9
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
tilman.metz's avatar
tilman.metz committed
11 12 13 14 15
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
Florian Wittkamp's avatar
Florian Wittkamp committed
16
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
tilman.metz's avatar
tilman.metz committed
17 18 19 20 21 22 23 24 25 26 27 28
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
*   Calculating source signal at different source positions with different
*   time-shift, centre frequency and amplitude (as specified in SOURCE_FILE).
*   Source signals are written to array signals 
*
*  ----------------------------------------------------------------------*/

#include "fd.h"


29
float ** wavelet(float ** srcpos_loc, int nsrc, int ishot, int SH, int STF){
30 31 32
    
    
    /* extern variables */
33
    extern int SOURCE_SHAPE,SOURCE_SHAPE_SH, NT, MYID,VERBOSE;
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
    extern float  DT;
    extern char SIGNAL_FILE[STRING_SIZE];
    extern char SIGNAL_FILE_SH[STRING_SIZE];
    extern FILE *FP;
    
    /*local variables */
    int nts, nt, k;
    float *psource=NULL, tshift, amp=0.0, a, fc, tau, t, ts, ag;
    float ** signals;
    
    
    
    /* ---------------------------- */
    /*   P SV Source                */
    /* ---------------------------- */
    if (SH==0) {
50 51
        if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
        if (SOURCE_SHAPE==7){
52
            psource=vector(1,NT);
53
            inseis_source_wavelet(psource,NT,ishot,SH,STF);}
54 55 56 57 58 59 60 61 62 63 64 65
        
        signals=fmatrix(1,nsrc,1,NT);
        
        for (nt=1;nt<=NT;nt++){
            t=(float)nt*DT;
            
            for (k=1;k<=nsrc;k++) {
                tshift=srcpos_loc[4][k];
                fc=srcpos_loc[5][k];
                a=srcpos_loc[6][k];
                ts=1.0/fc;
                
66
                switch (SOURCE_SHAPE){
67 68 69 70 71 72 73 74 75 76 77
                    case 1 :
                        /* New Ricker Wavelet, equal to SOFI2D */
                        tau=PI*(t-1.5*ts-tshift)/(ts);
                        amp=(((1.0-2.0*tau*tau)*exp(-tau*tau)));
                        break;
                    case 2 :
                        if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
                        else amp=((sin(2.0*PI*(t-tshift)*fc)
                                   -0.5*sin(4.0*PI*(t-tshift)*fc)));
                        break;
                    case 3 :
78
                        /* source wavelet from file SOURCE_FILE */
79 80
                        if (nt<=nts) amp=psource[nt];
                        else amp=0.0;
81
                        break;
82
                    case 4 :
83
                        /* sinus raised to the power of three */
84
                        if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
85
                        else amp=pow(sin(PI*(t-tshift)/ts),3.0);
86
                        break;
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
                        
                        break;
                    case 5 :
                        /* first derivative of a Gaussian */
                        ts=1.2/fc;
                        ag  = PI*PI*fc*fc;
                        amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts));
                        break;
                    case 6 :
                        /* Bandlimited Spike */
                        amp=0.0;
                        if(nt==1+iround(tshift/DT)){
                            amp = 1.0;}
                        break;
                    case 7 :
                        /* source wavelet from file SOURCE_FILE */
                        amp=psource[nt];
                        break;
                    case 8 :
                        /* integral of sinus raised to the power of three */
                        if (t<tshift) {
                            amp=0.0;}
                        if ((t>=tshift) && (t<=(tshift+ts))){
                            amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));}
                        if (t>(tshift+ts))
                        {amp=ts/(0.75*PI);}
                        break;
                    default :
115
                        declare_error("Which source-wavelet ? ");
116 117 118 119 120 121 122 123 124 125 126 127 128
                }
                
                
                signals[k][nt]=amp*a;
            }
        }
    
        if(VERBOSE){
            fprintf(FP," Message from function wavelet written by PE %d \n",MYID);
            fprintf(FP," %d source positions located in subdomain of PE %d \n",nsrc,MYID);
            fprintf(FP," have been assigned with a source signal. \n");
        }
        
129
        if (SOURCE_SHAPE==3 || SOURCE_SHAPE==7) free_vector(psource,1,NT);
130 131 132 133 134 135 136 137
        
        return signals;	
        
    }
    
    /* ---------------------------- */
    /*   SH Source                  */
    /* ---------------------------- */
tilman.metz's avatar
tilman.metz committed
138

139
    if (SH==1) {
140 141 142 143 144
        
        if (SOURCE_SHAPE_SH==3) {
            psource=rd_sour(&nts,fopen(SIGNAL_FILE_SH,"r"));
        }
        
145
        if (SOURCE_SHAPE_SH==7){
146
            psource=vector(1,NT);
Florian Wittkamp's avatar
Florian Wittkamp committed
147
            inseis_source_wavelet(psource,NT,ishot,SH,STF);
148
        }
149 150 151 152 153 154 155 156 157 158 159 160
        
        signals=fmatrix(1,nsrc,1,NT);
        
        for (nt=1;nt<=NT;nt++){
            t=(float)nt*DT;
            
            for (k=1;k<=nsrc;k++) {
                tshift=srcpos_loc[4][k];
                fc=srcpos_loc[5][k];
                a=srcpos_loc[6][k];
                ts=1.0/fc;
                
161
                switch (SOURCE_SHAPE_SH){
162 163 164 165 166 167 168 169 170 171 172
                    case 1 :
                        /* New Ricker Wavelet, equal to SOFI2D */
                        tau=PI*(t-1.5*ts-tshift)/(ts);
                        amp=(((1.0-2.0*tau*tau)*exp(-tau*tau)));
                        break;
                    case 2 :
                        if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
                        else amp=((sin(2.0*PI*(t-tshift)*fc)
                                   -0.5*sin(4.0*PI*(t-tshift)*fc)));
                        break;
                    case 3 :
173
                         /* source wavelet from file SOURCE_FILE */
174 175
                        if (nt<=nts) amp=psource[nt];
                        else amp=0.0;
176
                        break;
177
                    case 4 :
178
                        /* sinus raised to the power of three */
179
                        if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
180
                        else amp=pow(sin(PI*(t-tshift)/ts),3.0);
181
                        break;
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
                        
                        break;
                    case 5 :
                        /* first derivative of a Gaussian */
                        ts=1.2/fc;
                        ag  = PI*PI*fc*fc;
                        amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts));
                        break;
                    case 6 :
                        /* Bandlimited Spike */
                        amp=0.0;
                        if(nt==1+iround(tshift/DT)){
                            amp = 1.0;}
                        break;
                    case 7 :
                        /* source wavelet from file SOURCE_FILE */
                        amp=psource[nt];
                        break;
                    case 8 :
                        /* integral of sinus raised to the power of three */
                        if (t<tshift) {
                            amp=0.0;}
                        if ((t>=tshift) && (t<=(tshift+ts))){
                            amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));}
                        if (t>(tshift+ts))
                        {amp=ts/(0.75*PI);}
                        break;
                    default :
210
                        declare_error("Which source-wavelet ? ");
211 212 213 214 215 216 217 218 219 220 221 222
                }
                
                
                signals[k][nt]=amp*a;
            }
        }
        if(VERBOSE){
            fprintf(FP,"\n Message from function wavelet SH written by PE %d \n",MYID);
            fprintf(FP," %d source positions located in subdomain of PE %d \n",nsrc,MYID);
            fprintf(FP," have been assigned with a source signal. \n");
        }
        
223
        if (SOURCE_SHAPE_SH==3 || SOURCE_SHAPE_SH==7) free_vector(psource,1,NT);
224 225 226 227
        
        return signals;	
        
    }
228 229
    declare_error("Wrong input argument wavelet.c");
    return NULL;
tilman.metz's avatar
tilman.metz committed
230
}