• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Forward modelling of the Cotton-Mouton effect polarimetry on EAST tokamak

    2024-04-06 07:16:08MinyongSHEN沈敏勇JiboZHANG張際波YaoZHANG張耀YinxianJIE揭銀先HaiqingLIU劉海慶JinlinXIE謝錦林andWeixingDING丁衛(wèi)星
    Plasma Science and Technology 2024年3期
    關(guān)鍵詞:劉海衛(wèi)星

    Minyong SHEN (沈敏勇) ,Jibo ZHANG (張際波) ,Yao ZHANG (張耀) ,Yinxian JIE (揭銀先) ,Haiqing LIU (劉海慶) ,Jinlin XIE (謝錦林) and Weixing DING (丁衛(wèi)星),*

    1 School of Nuclear Science and Technology,University of Science and Technology of China,Hefei 230026,People’s Republic of China

    2 Institute of Plasma Physics,Chinese Academy of Sciences,Hefei 230031,People’s Republic of China

    Abstract Measurement of plasma electron density by far-infrared laser polarimetry has become a routine and indispensable tool in magnetic confinement fusion research.This article presents the design of a Cotton-Mouton polarimeter interferometer,which provides a reliable density measurement without fringe jumps.Cotton-Mouton effect on Experimental Advanced Superconducting Tokamak (EAST) is studied by Stokes equation with three parameters (s1,s2,s3).It demonstrates that under the condition of a small Cotton-Mouton effect,parameter s2contains information about Cotton-Mouton effect which is proportional to the line-integrated density.For a typical EAST plasma,the magnitude of Cotton-Mouton effects is less than 2 π for laser wavelength of 432 μm.Refractive effect due to density gradient is calculated to be negligible.Time modulation of Stokes parameters (s2,s3) provides heterodyne measurement.Due to the instabilities arising from laser oscillation and beam refraction in plasmas,it is necessary for the system to be insensitive to variations in the amplitude of the detection signal.Furthermore,it is shown that non-equal amplitude of X-mode and O-mode within a certain range only affects the DC offset of Stokes parameters (s2,s3) but does not greatly influence the phase measurements of Cotton-Mouton effects.

    Keywords: EAST,Cotton-Mouton effect,polarimeter interferometer,electron density measurement,Stokes vector

    1.Introduction

    The purpose of this paper is to propose a method for continuous and reliable high-precision measurements of plasma electron density [1] using a polarimeter interferometer [2]based on the Cotton-Mouton effect (hereafter abbreviated as CM) on EAST tokamak or on future fusion devices.This measurement can provide key information for density control of magnetic fusion reactors [1].Unlike conventional interferometers that require time accumulation for phase wrap-up,the CM effect induces a phase shift of approximately 15 degrees on the central chord of EAST when a wavelength of 432 μm is chosen,effectively compensating for the disadvantage of fringe jumps when the phase difference exceeds 2 π after the interruption and reconnection of the interferometer signal [3].To guarantee the density accurate measurement,CM effect-based polarimetry can be operated simultaneously with an interferometer that employs the same probing beam [4].Furthermore,due to the instabilities arising from laser oscillation and beam refraction in plasmas,the system is designed to be insensitive to variations in the amplitude of the detection signal [5].

    There are two types of polarimeters: those based on the Faraday rotation (hereafter referred to as FR) effect and those based on the CM effect [6].These magneto-optical effects depend on the magnetic field aligned parallel or perpendicular to the propagation direction of electromagnetic waves in the plasma [4],respectively.Therefore,the choice of one of these effects depends on the optical beam path.For density measurement,the Faraday rotation effect is not suitable due to variation of magnetic field parallel to the beam path in poloidal section.Hence,the CM polarimeter becomes a good option for high density measurement without fringe jump [7].This paper limits the configuration to within a poloidal section.

    In the context of superimposed elliptical modulation,achieving exacting power equality between two orthogonal laser beams can be nearly impossible [8].This study aims to demonstrate that,even under this condition,accurate measurements of the relatively small CM effects can still be obtained through a phase comparator,while remaining insensitive to changes in the detection signal amplitude.Furthermore,we elucidate the rationale behind selecting radiation with wavelengths falling within the far-infrared(FIR) range.Although a Cotton-Mouton polarimeter has recently been developed on the Compact Helical System(CHS) [4] and is proposed for the International Thermonuclear Experiment Reactor (ITER) [9],its applications are fewer than those of Faraday polarimeters.We will develop a CM effect polarimeter system that utilizes a 432.6 μm formic acid laser on EAST tokamak.In section 2 we introduce typical parameters of EAST plasma and present some optical quantities related to polarimetry.In section 3,we utilize Stokes equation and Mueller matrix to investigate evolution of polarization of lasers in EAST plasma.By using time modulation of Stokes parameters,one can determine elements of Mueller matrix which relates to plasma parameters.In section 4,we study the coupling between FR effect and CM effect.An explicit nonlinear relationship is proposed between CM phase shift and χ for phase difference not much less than one.Section 5 provides a brief summary.

    2.The propagation of electromagnetic waves on EAST tokamak

    2.1.EAST plasma

    EAST is a magnetic confinement tokamak fusion device,with a plasma major radius of approximately 1.9 m,an average minor radius of approximately 0.45 m,and a magnetic field strength of about 2.4 T along the magnetic axis.Currently,there are five vertical ports in EAST device,evenly spaced at intervals of 90 mm in theR-direction,located atR=1640 mm,1730 mm,1820 mm,1910 mm,and 2000 mm,respectively.The vertical port atR=1910 mm approximately passes through the magnetic axis.The laser beam is injected from the upper port and detects the lineaveraged electron density at the lower window of the device.

    The magnetic configuration of EAST shot #98958 at 21 s is shown in figure 1,where the red dashed line represents the last closed magnetic surface.A right-handed orthogonal coordinate system is established,with the vertical direction defined as theZ-axis,the direction from the high-field side to the low-field side in the major radial direction as theXaxis,and the toroidal direction as theY-axis.In the simulations performed in this study,the incident position is set within the range ofR0∈(1640 mm,2000 mm).

    x0is defined as

    The vertical channel passing through the magnetic axis is defined asx0=0,which corresponds to the range ofx0∈(-270 mm,90 mm)in the minor radial direction.Electron density (ne) and toroidal magnetic field distribution(Bt) for EAST shots # 71320 at 3.3 s and #98958 at 21 s are shown in figures 2 and 3,respectively.

    2.2.Elliptical polarization of laser beam [8]

    In order to study the evolution of the polarization state of electromagnetic waves in magnetized plasmas,the following conventional quantities are introduced to describe the elliptical polarization state of a laser beam.

    A top view of elliptical polarization of laser beam is shown in figure 4,where beam propagates along the negative direction of theZ-axis defined in figure 1.

    In the electric field coordinate systemX-Y,α is the arctangent of the ratio of electric field amplitude in theYdirection to that in theXdirection.

    Generally,along the propagation direction,the projection trajectory of the polarization electric field vector’s endpoint onto theX-Yplane is an ellipse.Polarized ellipse in the elliptical coordinate systemX′-Y′,χ is the ellipticity,referred to as the arctangent of the ratio of the short axis of the ellipse to the long axis

    ψis the azimuthal angle of the polarization plane,referred to as the angle between the long axis of the ellipseX′and theX-axis in the electric field coordinate system.ψFRis the change in the angle of the polarization plane due to Faraday effect

    Plasma can be viewed as an anisotropic birefringent medium.The mode of the electromagnetic wave in which the electric field is parallel to the magnetic field direction is the O-mode,while the X-mode is perpendicular to the magnetic field.The passage of linearly polarized light through birefringent medium results in ellipsoidalization,known as the Cotton-Mouton effect.

    Upon entering a plasma,incident polarization states will all split into fast and slow modes,each with distinct optical path lengths upon exiting.The term “fast” indicates a smaller optical path length.The slow state is delayed relative to the fast state,and this delay is due to the difference in optical path lengths,which describes the relative phase change between the two states [10].

    As shown in figure 5,a retarder splits light into two linearly orthogonal polarized components.Half-wave linear retarder delays one of the linearly polarized components by half a wavelength and is useful for changing the direction of linearly polarized light.With a vertical fast axis (lower refractive index),the component perpendicular to the magnetic field direction is referred to as the ‘X-mode’.The X-mode enters the half-wave plate along the fast axis and undergoes 2 and 1/4 cycles.The O-mode,on the other hand,traverses 2 and 3/4 cycles along the slow axis,resulting in a phase delay of π.Their combined light transforms from rightcircularly polarized light to left-circularly polarized light.

    As shown in figure 4,an inscribed circle can be constructed using the length of the minor axis of the ellipse as its diameter.Constructing a line parallel to lineX′and passing through point A,let AC be parallel toX′.Similarly,constructing another line parallel to lineX′and passing through point B,let BD be parallel toX′.These lines intersect the inner circle at points C and D,respectively.The angle ∠COD is geometrically referred to as the eccentric anomaly and is denoted as δ.It is also known as the orthogonal electric vector phase. δCMis the change in the orthogonal electric vector phase angle due to the CM effect [11].

    where the constantCCMB⊥is perpendicular to the direction of light propagation of the magnetic field.WhenB⊥,ne,λ andlare expressed in SI,CCMhas a numerical value of 2.4568×10-11,the unit of δCMis rad.

    Figure 6 illustrates the phase shift of the CM effect in theX-direction for each chord.The maximum phase shift of the CM effect does not exceed 20°,thus avoiding the occurrence of fringe jumps that would result from exceeding 2 π.

    To achieve a phase shift magnitude of approximately 20°the wavelength of the far-infrared laser was chosen as 432.6 μm.According to equation (5),figure 7 illustrates that the phase shift is proportional to the cube of the wavelength.Due to the refractive effect,the plasma region traversed by the same chord slightly varies at different wavelengths,so the curve in this figure is not strictly a cubic function.

    Plasma can exhibit refractive effects due to density gradients.The propagation of light,as described by Born in[12],is governed by Snell’s Law,which relates the incident angle (θ1) and the refracted angle (θ2) to the absolute refractive indices of the two media,denoted asn1andn2:

    This law describes the relationship between the angles of incidence and refraction when light passes from one medium to another.By calculating the refractive index vector at each spatial step and determining the displacement after passing through the plasma,the quasi-optical deviation caused by the presence of plasma is inferred.This displacement must be kept as small as possible so that probing beam can exit from vacuum windows for signal detection.

    119 μm represents the CH3OH laser,195 μm represents the DCN laser,337 μm represents the HCN laser,432.6μm represents the CH3COOH laser,and 462 μm represents the solid source.figures 8 and 9 use data from shot #71320,which can be observed that the wavelength does not change the trend of lateral offset and phase shift.However,due to the proportional relationship between phase shift and λ3,the wavelength significantly affects the magnitude of the measured values.When the wavelength is smaller than the chosen wavelength of the formic acid laser in this paper,the peak phase shift is less than 10°,which will affect the resolution of the polarimeter interferometer [13].

    Figure 9 illustrates the variation in the magnitude of transverse displacement caused by the refractive effects in each channel along theX-direction,which is on the order of millimeters.After consulting the engineering drawings,it was determined that the diagnostic window has a diameter of 70 mm.During the discharge,the laser generates a refraction of approximately 2 mm,and the far-infrared laser spot diameter is around 30 mm.As long as the optical alignment is accurate on the tabletop during experiments,it can be observed that the currently selected wavelength meets the requirements of the diagnostic window size in EAST,enabling the reception of laser signals emitted in the opposite direction.

    When the probe beam is incident vertically,B⊥is equivalent toBt,which is approximately constant along the∫ light path,as shown in figure 3.The density chord integraldlis measured using the interferometer at EAST.The CM effect along this chord is calculated by transforming it through equation (5),as δCM=so the measurement of electron density along the central chord of the poloidal section in the plasma is most appropriately performed using a vertical port CM effect polarimeter interferometer.It can be seen that the phase shift of the CM effect is expected to be around 15° throughout the EAST discharge,a vertical-view Cotton-Mouton polarimeter interferometer is capable of providing reliable density measurement.In figure 10,there is a sharp and brief density jump att=73 s,which is caused by the interferometer’s fringe jump,and the CM effect polarimeter interferometer will effectively circumvent this anomaly.

    The CM effect acts directly on the phase shift change dδCM,which leads to the ellipticity χ change through equation (7)[8].

    Figure 11 illustrates the noticeable variations in ellipticity along the central chord,demonstrating the theoretical feasibility of measuring the CM effect through the ellipticity information carried by the polarized ellipse.

    Figure 1.Diagram of the definition of the coordinate system.Blue line represents the flux surface.Red dashed line represents the last closed magnetic surface.

    Figure 2.The typical distribution of ne in the poloidal section for EAST shots #71320 and #98958,where ρ represents the normalized magnetic flux.

    Figure 3.The typical distribution of Bt in the poloidal section for EAST shots (a) # 71320 and (b) #98958.Display only the portion of normalized magnetic flux that is less than or equal to 1.The red arrow indicates the approximate location of viewing chords.

    Figure 4.Schematic representation of α,δ,χ and ψ in the polarization ellipse.

    Figure 5.A schematic representation of polarized light entering a birefringent medium -a linear half-wave plate,depicting the phase delay.Three states of polarized light from top to bottom are the the O-mode,X-mode,and their combined light.

    Figure 6.Phase shift of CM effect in the X-direction channels,where the horizontal axis represents the position of the incident laser.The solid line represents shot #98958,and the dashed line represents shot #71320.

    Figure 7.The relationship between the phase shift of the CM effect along the central chord and the wavelength of far-infrared laser.The black dot represents the selected wavelength of 432.6 μm.

    Figure 8.The relationship between phase shift and incident position under different laser wavelengths (shot #71320).

    Figure 9.The relationship between the incident position and the lateral offset under different laser wavelengths (shot #71320).

    Figure 10.CM effect transformed with time,where the string integral of the density is given by a solid source interferometer arranged in the centre channel.The Bt approximation is considered equal along the course and is regarded as the magnetic field at the magnetic axis.The wavelength is taken to be 432.6 μm (shot 98958)

    Figure 11.The schematic diagram depicting the variation of the polarized ellipse along the central chord,where the abscissa represents Ex and the ordinate represents Ey.

    Figure 12.Stokes parameters map the polarization state to a point on the Poincaré sphere,which is a unit sphere in an abstract 3-D space.

    Two concepts should be clarified here.The term “orthogonal” in the context of orthogonal electric vector phase refers to the perpendicular relationship in theX-Ycoordinate system of the electric field.The polarization ellipse is the trajectory traced by the endpoints of the electric field vector after one cycle,i.e.,a 360° change in the phase angle.From the extremal point B ofExto the extremal point A ofEythrough a certain phase angle,which is the mapped ∠COD,not necessarily 90°,so-called “quadrature”.CM effect δCMacts directly on the change of phase angle of the orthogonal electric vector,which is a concept of polarization optics.The“phase angle” here is not equivalent to the phase change of the light intensity.

    3.Evolution of polarization state based on Stokes equation

    As discussed above the measurement of polarized ellipse is non-trivial,it requires two measurements ofExandEy.Accurately measuring the absolute electric field amplitude can be very challenging.Instead the signal is modulated to incorporate the information of the CM effect within the phase difference of the Stokes parameters.One can modulate the polarized ellipse in time so that phase difference betweenExandEycan be determined.In this section we utilize Stokes equation,Jones vector and Mueller matrix to describe evolution of polarization state in EAST.By detecting Stokes parameters one can determine plasma density,magnetic field.

    3.1.Description of the polarization state of electromagnetic waves based on the Poincaré sphere

    The purpose of introducing Stokes parameters is to provide a mathematical representation of the polarization state of electromagnetic waves.Stokes parameters encompass polarization characteristics such as intensity,degree of polarization,and polarization ellipse orientation.They enable a quantitative description and analysis of the polarization state of light[14].

    Introducing the Poincaré sphere with a unit radius of 1 in the (s1,s2,s3) space facilitates the study of polarization,as depicted in figure 12.Here,all possible polarization states(ψ,χ) are uniquely mapped by points P on the unit sphere surface.Establishing an orthogonal Cartesian coordinate system with the center of the sphere as the origin,the three directional components can be represented using a one-rowthree-column matrix of Stokes parameters,corresponding to a latitude of 2 χ and a longitude of 2 ψ.

    The normalized Stokes parameters can be defined as [15]

    Two orthogonal polarizations are represented by antipodal points on the equator of the Poincaré sphere.Taking the X and O beams under investigation in this paper as examples,the horizontal polarization is characterized by the Stokes parameterssH=(1,0,0),and the vertical polarization is characterized bysV=(-1,0,0).When polarizations propagate through a homogeneous non-absorbing medium,they experience different characteristic refractive indices,denoted as μ1and μ2,respectively,satisfying μ1>μ2.

    Ωrepresents the angular velocity of the rotation of the Stokes parameters on the Poincaré sphere [16].The fast axis aligns with the direction of the Ω vector,passing through the two points on the Poincaré sphere that represent the characteristic orthogonal polarizations,sHandsV.An introduction to the fast axis can be found in figure 5.Therefore,the evolution of polarization can be expressed by a vector equation [17].

    Here,the magnitudes ofCFRandCCMare given by equation(4) and equation (5),respectively; Ω1,2,3represent the angular velocities of the Stokes parameters rotating around thes1,2,3axes,and each component is contributed by the CM effect dδCM,cross-terms,and the FR effect dψFR.The quantities Ω1,2,3reflect the optical properties of the plasma medium.By setting an appropriate incident laser polarization direction and measuring the changes in polarization optical quantities,Ω1,2,3can be determined,providing information about plasma density,magnetic field,and other parameters.

    Figure 13 illustrates when the central chord of the light beam varies along the negative direction of theZ-axis,the impact of the CM effect on three components of the Stokes parameters.The initial Stokes parameters are given by the superposition of horizontally and vertically linearly polarized light,i.e.,s(z=0)=(0,1,0).The variation in the amplitude ofs3is -0.2444 rad,which can be determined through the ∫ vector cross product in equation (10) . Since,the CM effect on EAST can be approximated ass3=0.2508 rad,relating plasma parameter.

    Figure 13.The variation of the Stokes parameters along the central chord when the central chord of the light beam varies along the negative direction of the Z-axis.

    Figure 14.Addition of equal horizontally and vertically polarized beams with a 20 % frequency difference yields a time-varying polarization state [10].

    Figure 15.Schematic diagram of elliptical polarization modulation trajectory under different transmission ratio conditions by superimposing X-mode and O-mode.

    Figure 16.Comparative analysis of the time modulation of s2 (equation (27)) and the simulation results based on definition (equation(10)).The comparative analysis of the simulated results for time modulation of s2.The left plot represents equal power ratio transmission,while the right plot represents non-equal power transmission.

    Figure 17.In the case of a significant CM effect (×3),comparative analysis of the time modulation of s2 (equation (27)) and the simulation results based on definition (equation (10)).The left plot represents equal power ratio transmission,while the right plot represents nonequal power transmission.

    Figure 18.Schematics of the optical arrangement of the CM polarimeter.

    3.2.Derivation of elliptical modulation based on Jones vectors

    For application of polarized ellipse time modulation technique,we present the time-dependent Stokes parameter using Jones vectors.From equation (8),it can be observed that the Stokes parameters represent the temporal ensemble average of the orthogonal components of the electric field and cannot reflect the frequency modulation of a single beam of light.Therefore,the Jones vector representation is introduced,which employs a 2 ×1 matrix to represent thexandycomponents of the electric field [17],

    The electric field resulting from the superposition of the Xmode and O-mode is

    where,Δω=ω2-ω1represents the intermediate frequency of the X-mode and O-mode,and the subscript ‘a(chǎn)’ denotes the amplitude,and δ1and δ2are initial phases for each component.From equation (2),the amplitude ratio α can be obtained as

    The parameter α is a constant that is solely related to the ratio of the transmitted power. Furthermore,sinces1=cos(2α),the trajectory of elliptical modulation is reflected on the Poincaré sphere as a constant value ofs1.Phase difference

    the trajectory of the modulation process on the Poincaré sphere is represented by a circle drawn from the periodic variation of the phase angle.

    In figure 14,the end of the electric field vector traces a more complex shape,which can be interpreted as a timevarying polarization state resembling an evolving ellipse.This distinctive shape is known as a Lissajous figure,which is formed by the parametric equation.To a Stokes polarimeter,this beam appears as unpolarized light and cannot be distinguished.

    From equation (8),the Stokes parameters of the synthesized light resulting from the superposition of horizontal and vertical polarization modulation can be obtained as

    Non-uniform power emission will introduce deviations in the intensity amplitude of the measurement by affecting the magnitude ofs1.The elliptical modulation trajectory under different transmission ratio conditions is shown in figure 15,represented by a series of circles with their centers falling on thes1axis.The standard form in an orthogonal Cartesian coordinate system between the variabless2ands3can be written as:

    which is represented by the yellow dashed line.In general,there exists a nonlinear relationship between Δχ and δCM.

    Specifically,when two orthogonally polarized beams are equally transmitted with equal power,i.e.,EaX=EaOor tanα=1,from equation (7) it can be deduced that δ and ellipticity χ have a linear relationship.The relationship between the variabless2ands3is defined by the equationsin2(2×45°),forming a meridian loop byψ=90°and ψ=-90°.The circle can be analogously compared to geographical east longitude 90° and west longitude 90°,as i ndicated by the purple dashed line in figure 15.

    3.3.Approximate analytic solutions ofs2

    The output polarization state can be characterized using either Stokes parameters or,equivalently,ψ and χ.The former is often more convenient since the detector signals are directly related to the three components of the Stokes parameters.For the design and analysis of polarization measurements,the Mueller matrixMis required to comprehensively describe the effects of optical components or plasma media on the polarization state [18].This matrix is the characteristic of plasma (density and magnetic field etc).The output Stokes vector relates to Mueller matrix (plasma parameters) as the following expression [17]:

    Due to the non-commutativity of the Mueller matrices for the two effects,it is not possible to simply integrate each effect separately to calculate the total effect.Therefore,when both effects are strong,the coupling cannot be separated[19].When at least one of the two effects is small,approximate analytic solutions of the evolution equation can be obtained.When the FR effect is small,we consider two cases for discussion: the CM effect where δCM?1 or δCM~1 .The dimensionless parametersW1,2,3are used for analysis [17].

    The parametersW1,2,3correspond to the CM effect δCM,the cross term,and the FR effect ψFRrespectively.(1)W1,2,3?1

    According to equation (18),Asatisfies the following equation:

    where now

    where the components of the matricesM(1),M(2),M(3)··· are given by

    where again repeated indices imply summation [20].

    When the plasma effect is small,the Mueller matrix can be expanded to second order,

    For a tokamak having up/down symmetry with respect to the equatorial plane,the following relations hold:W2=0,W13=W31=W12=-W21,andW23=-W32.These relations are investigated under the assumption that?1 and?1,therefore [17],

    Based on equations (16) and (18),and using the the auxiliary angle formula,the second components2of the Stokes parameters can be expressed as follows:

    So,multiplyingM21by the constant cos2α only affects the DC offset ofs2,while multiplyingM22andM23by terms involving (Δωt) affects both the amplitude and phase ofs2.

    The Stokes parameters2of laser polarization state after plasma is modulated at frequency (Δω).It contains information about the chord integralW1,The phase ofs2does not include quantities related to the power ratio α.One can determineW1by phase measurement instead of amplitude measurement (sin2α).In the process of orthogonal superposition signal radiation,the (Δωt) changes while α remains constant.The simplest case isEax=Eay,α=π/4,where two laser radiations have equal power.

    The errors resulting from non-equal power radiation is discussed here.Furthermore,parameters in Mueller matrix may not be very small.These approximations can introduce finite corrections.A comparison is made between the definitional formula (equation (10)) and the computational formula(equation (27)) fors2,and the results are plotted in figure 16.Without loss of generality,we set the modulation phase(Δωt)to be 90°,the phase variation of the AC component is reflected in a larger amplitude variation.

    In the case of non-uniform power emission,without loss of generality,the power ratio is assumed as α=Under equal power emission,the result from the definitional formula iss2d=0.252,while the result from the computational formula iss2e=0.248.Under non-uniform power emission,the result from the definitional formula iss2d=0.263,while the result from the computational formula iss2e=0.260.The errors introduced by neglecting higherorder terms in the computational formula are only -1.59%and -1.14% respectively.Now let us perform an order of magnitude analysis on the sources of error.As shown by equation (26),s2(z) consists of three terms,each containingM21,M22,andM23.The errors arise from neglecting terms of order three and higher,beyondM(3),in the derivation from equation (20) to equation (24).In the previous section,we assumed Δωt=90°,making the second term in equation (26)equal to 0.From equation (22) and equation (23),we can derive the expressions forwith respect toW

    In equation (26),sinceWis always positive,whenα increases,the weight ofM23becomes greater,leading to an overestimation ofs2by neglectingConversely,when α decreases,the weight ofM21increases,resulting in an underestimation ofs2by neglecting

    As previously mentioned,in the case of vertically incident CM effects on the poloidal section,the subscripts indicate that index 3 is attributed to CM effects and index 1 is attributed to FR effects. Therefore,|-W333|>|W111|,|-W322|>|W122|,|-W311|>|W133|.In other words,

    This indicates that the error introduced when α takes a smaller value must not be greater than that introduced when α=π/4.

    This indicates the feasibility of the proposed heterodyne measurement scheme with small errors.It can be understood that the non-equal power component only generates a DC offset of (W3-W21)cos2α and a signal attenuation of sin2α times at the mixer stage.It is important to emphasize that the derivation in this subsection cannot be extended to cases whereW1?1 is not satisfied.

    (2)W≈W1;W2,3?1

    For example,in the case of vertical propagation in a future tokamak,it may not always satisfy the condition|W1|?1in practical scenarios.Now let us write the analytical expression for the first-order expansion ofs2(z) whenW~O(1).

    Similarly to equation (26),the second components2of the Stokes parameters is given by [17]

    where,

    The derivation process of expandings2(z) to the first order whenW~O(1) is provided in Appendix A.

    To simulate a larger CM effect,the typical toroidal magnetic field magnitudeBtis multiplied by a factor ofeffectively amplifying the CM effect by a factor of 3.figure 17 is plotted to illustrate the measurement under this scenario.Without loss of generality,we set the modulation phase (Δωt) to be 90°,the phase variation of the AC component is reflected in a larger amplitude variation.The power ratio α in the case of unequal power emission isthe same as in figure 16.

    In the case of equal power emission,the result from the definition formula iss2d=0.694,while the result from the calculation formula iss2e=0.746.Similarly,in the case of unequal power emission,the definition formula yieldss2d=0.646,and the calculation formula yieldss2e=0.692.The errors introduced by neglecting higher-order terms in the calculation formula are +7.49 % and +7.12 % respectively.

    If a longer wavelength laser is chosen to achieve a larger CM effect,this heterodyne measurement scheme will introduce larger errors.Moreover,the longer the wavelength,the more pronounced the refractive effect becomes,so the detector no longer receives the laser signal precisely at the center,leading to a decrease in signal-to-noise ratio.Therefore,in practical experiments,the CM effect is typically not allowed to become excessively large (around δCM~O(1)) during electron density measurements.

    In summary,under the condition of non-equal power emission of orthogonal light,the superimposed elliptical modulation can still accurately measure relatively small CM effects through a phase comparator,while being insensitive to variations in the detected signal’s amplitude.In order to achieve an appropriate magnitude of the CM effect,the wavelength of the radiation used should fall within the far infrared (FIR) range.The error arises from neglectingand,so non-uniform power emission only affects the signal-to-noise ratio,with minimal impact on the calculation error of the CM effect phase.

    3.4.Determination of Stokes vector

    Elliptical modulation enables the direct measurement of δCMthrough the phase difference of the Stokes parameters.In this section,we will discuss the polarization optical components and their configurations in the polarization analyzer system.From equation (16),boths2ands3contain terms modulated with respect to (Δωt).

    The intensity measurement ofs2can be performed without the presence of a waveplate.In principle,it is feasible to measures3as well.However,fors3measurement,a quarterwave plate needs to be installed in the optical path,requiring prior measurement of the waveplate’s absorption factor which involving additional effort [21].Therefore,we ultimately decide to measures2.

    As shown in figure 18,the micro-displacement platform is used to adjust the length of the resonant cavity,ensuring that the frequency difference between the two laser sources is equal to the desired beat frequency.Modulating the orthogonal light generates an intermediate frequency signalIF~O(1 MHz),facilitating heterodyne measurement.One laser source emits linearly polarized light,which is passed through a half-wave plate to be orthogonal to the polarized light emitted by the other laser source.The beams are then combined,adjusted to a collinear state,and guided through an optical waveguide to the top of EAST device,where they are incident on the plasma.

    Direct measurement cannot determine the value of a specific component.The power measured by the detectorID,is the weighted sum of the individual components of the Stokes parameters [8].

    where,I0represents the intensity before passing through the polarizer,θ is the angle between the transmission axis of the polarizer and theX-axis.We set θ=45° to ensure that the intensity only contains information abouts2,

    Under the assumption of a constantI0,there is a linear relationship between the intensity ands2,withs2having the highest coefficient,thus the phase measurement ofs2can be transformed into a phase measurement of the intensity.Due to the difference in refractive index between the two beams in the plasma,the phase shift difference between the ordinary wave and the extraordinary wave is compared using a mixer and a phase discriminator.This allows for the determination of the chord-integrated plasma chord density.

    4.Measurement of Cotton-Mouton effect based on ellipticityχ

    When the phase difference of the CM effect is only about 15 degrees,reducing phase noise becomes an important issue.Additionally,when the CM effect is significant,heterodyne measurement methods can introduce an approximate 10%measurement error.Now,we will demonstrate that an accurate analytical solution for the ellipticity χ is feasible under larger Cotton-Mouton (CM) effects.

    4.1.The system of differential equations involving ψ andχ

    Stokes parameters are observable,but they are not as intuitive as the polarization ellipse parameters.To gain a more comprehensive understanding of the changes in the polarization ellipse within the plasma,it is advantageous to convert the Stokes equations into expressions based on the parameters of the polarization ellipse.

    To express the differential interaction between χ andψ[22]

    During the propagation of light in a plasma,the presence of the Faraday rotation (FR) effect causes the azimuthal angle to vary continuously.Consequently,sin2ψ is not a constant.By equation (35),ψ is coupled explicitly to χ.The FR effect,although not acting directly on χ,couples implicitly to χ by changing ψ through equation (36) [23].

    4.2.The expression ofψ(χ)

    If one of the FR and CM effects can be neglected,the other effect,no matter how strong,can be written as an explicit expression and can theoretically be measured [20].When the CM effect can be ignored,the FR effect can be measured no matter how strong it is,and this conclusion was verified in JET [24].This conclusion is obviously.When δCM=0,equations (35) and (36) are written as

    Conversely,when the FR effect can be neglected and the CM effect plays a dominant role,the problem can be considered as a generalized Volterra-like problem with nonconstant coefficients [25].

    equation (40) is divided by equation (39) to obtain

    Because the indefinite integral

    wherexis an arbitrary function andCis a constant.

    Therefore,integrating equation (42) to obtain,

    where χ0and ψ0are known constant values representing the polarization parameters of the incident light.

    When |cos2ψ0|·cos2χ0≠0,it satisfies

    Furthermore,define Const as

    According to the properties of trigonometric and inverse trigonometric functions,the constant Const ∈(0,1] in order to ensure that the value domain of ψ remains within the set of real numbers.

    Figure 19 depicts the relationship between ψ andχ obtained from equation (49),where the arrows indicate the direction of phase space trajectories influenced by the CM effect.Points p1and p3represent special cases at χ=0,where ψ reaches extremum.Similarly,points p2and p4correspond to ψ=0,where χ attains extremum.

    Figure 19.Considering different values of Const,the phase space trajectories of χ-ψ and the extremal points p1-4 are depicted.The arrows indicate the direction influenced by the CM effect.

    Figure 20.The integrand of dδCM (equations (54) and 63) under different Const,where p1-4 is the same as defined in figure 19.The abscissa of the black vertical dashed line represents the inflection point χ=A.Different line types correspond to different emission power ratios.The blue solid line represents equal power emission with Const=0;the green dashed line represents Const=;and the red dotted line represents Const=.The case of Const =1,which corresponds to point A coinciding with χ=0,cannot be depicted in this figure.

    4.3.The explicit expression after decouplingδ(χ)

    If the initial polarization parameters (χ0,ψ0) are given and the end-state χ is measured,the phase shift of the CM effect can be written as a function dδCM=dδCM(χ0,ψ0,χ),which achieves the decoupling of χ from ψ.

    Noting that at the singularity where ψ=0,the denominator on the right-hand side of equation (41) becomes 0 in the improper integral.EAST satisfies thatBzis homogeneous along the range in thezdirection and ψ varies monotonically along the range.Without loss of generality,ψ is assumed to be monotonically decreasing.We classify and discuss two distinct scenarios for the phase space trajectories based on the positive and negative values of ψ.The first scenario involves trajectories moving from p4to p1and then to p2,while the second scenario entails trajectories moving from p2to p3and then to p4.

    (1) ψ in the range greater than 0.Equation (47) is written as

    substituting equation (35) into equation (50),

    According to equation (51),the integrand function can be computed as follows:

    The definite integral is

    The integral calculation procedure of equation (55) is shown in Appendix B.

    (2) ψ in the range less than 0.

    If the polarized light propagates such that ψ=0,i.e.,(cos22χ-Const)=0,the singularity of this elliptic degreeχ can be defined as

    by equation (42),

    from equation (35),when ψ<0,then sin2ψ<0.Implying d χ<0,χ is monotonically decreasing,i.e.,χ

    holds true,whereCis an arbitrary constant.Introducing the terms=A-χ > 0 as the difference between the final state of ellipticity and the value ofA,we can express d χ=-ds.Substituting equation (58) into equation (57) results in:

    The integrand function can be computed as follows:

    According to equations (54) and (63),the phase space trajectories of the integral function family for different values of Const are illustrated in figure 20,where the ordinate represents the dimensionless quantity “Integrand”.The vertical dashed line represents the inflection point at ψ=0,denoted as χ=A.Following the assumptions stated earlier,when χ>0,ψ monotonically decreases,resulting in clockwise phase space trajectories.

    Based on the aforementioned analysis,the following conclusions can be inferred:

    (1) The area enclosed by the “Integrand” curve and theY-axis (y=0) represents the definite integral of the CM effect cumulating along the trajectory.This area has a maximum range of π and is computed as the definite integral with the limits of χ ranging from -AtoA.

    (2) We will now explore certain scenarios with distinct initial values of Const.

    (i) In the process of solving the system of differential equations,Const ≠0 was employed.Specifically,the values of (χ0,ψ0) are set such that Const=0,resulting in a simplification of equation (55) as follows:

    Based on figure 20,it is evident that in this particular case,the domain of χ exhibits the widest range.

    (ii) In the scenario where Const=1,the domain of (χ0,ψ0) is restricted to a single isolated point at (0,0),indicating that the incident laser is solely a horizontally polarized light.This isolated point,illustrated in figure 19,corresponds to point A coinciding with χ=0,which cannot be depicted in figure 20.Consequently,this confirms the inference that measuring X-mode or O-mode alone is insufficient.

    (3) When Const ≠1,the CM effect can be explicitly expressed in terms of the initial and final states of elliptical parameters.Notably,assuming the FR effect can be disregarded,this decoupling occurs independently of the final state solution for ψ.

    4.4.Commonalities with heterodyne measurements

    Now,we consider the relationship between non-uniform power emission superposition modulation and ellipticity measurement from the perspective of Poincaré sphere.The simultaneous combination of the definition equation (48) for Const and the definition equation (8) for Stokes parameters yields,

    This indicates that both homodyne measurement and heterodyne measurement of orthogonal light with arbitrary power ratios share a common characteristic ofs1(0) remaining constant over time.This can be derived from equation (9) [20].

    When the modulation process crosses the position ofs2=0,the monotonicity of χ changes,which corresponds to the singularities mentioned in section 4.3.If=0,the modulation process is represented by the purple dashed line in figure 15,and the value range of the measured parameters2is maximized.This implies the highest resolution for equalpower transmission of orthogonal light.If=1,the modulation process degenerates into an isolated point H as shown in figure 15,wheres2≡0,meaning that there is no vertically polarized light emitted,making it impossible to measure the CM effect.

    4.5.Modeling results and discussion of ellipticity measurement

    In order to simulate phase shift measurements under higher parameter conditions,a larger Cotton-Mouton (CM) effect than that of existing tokamak devices is considered.By comparing the definition equation (5),we verify the accuracy of the nonlinear calculation formula equation (55) for the CM phase shift δCMin ellipticity measurements and demonstrate its superiority over the previous linear approximation formula equation (65).

    Figure 21 presents the results for the same shot (shot#98958) at the same moment (21 s),where the toroidal magnetic fieldBtis multiplied byeffectively amplifying the CM effect by a factor of 3.The vertical axis represents the phase shift,and the horizontal axis represents the incident positionx0.Different incident elliptical polarization conditions,referred to as State 1,State 2,and State 3,are considered to simulate a more general initials1.State 1 corresponds to the initial polarization parameters (χ0,ψ0) ofState 2 corresponds toand State 3 corresponds toIn State 1,State 2,and State 3,the values ofs1are 0,0.309,and 0.294,respectively.

    Figure 21.Comparison of CM phase shift calculation results between nonlinear computational formula and the definition formula under different incident elliptical polarization conditions.The results calculated using equation (5) are represented by dashed lines,the results calculated using the nonlinear calculation formula equation (55) are represented by “ ×”,and the results calculated using the previous linear approximation formula equation (65) are represented by “+”.

    The error is defined as

    where ? represents the phase shift.

    When the polarization parameters (χ0,ψ0) arewhich corresponds tos1=0,the error of the approximation formula is almost zero,whereas for other cases,the error becomes more pronounced.The error reaches 7.32 % for State 2 and 11.55 % for State 3.The computational equation is perfectly fit with the phase shift calculation of the defined equation in either incidence state,and the error in the numerical calculation is small,about 0.5 %,which is basically negligible.

    In summary,the three subplots in figure 21 correspond to the scenarios in equation (18),where the plasma Mueller matrixM(z) is the same,but the initial Stokes parameters(0)is different.Since the plasma state remains the same throughout the passage,if the nonlinear computational model given by equation (55) is correct,the resulting δCMshould be consistent with the computed results obtained from equation (5).The nonlinear equation makes polarization measurements without having to reduce the wavelength and sacrifice resolution for the applicability of the linear model,ensuring the measuring accuracy of large CM effect.

    Next,we will discuss the second case mentioned in section 4.3,which corresponds to the situation when ψ crosses the zero point,i.e.,across p2and p4(as shown in figure 20).In this case,the change in χ does not exhibit monotonicity,resulting in the integrand function of dδCMfirst being greater than 0 and then becoming less than 0 after crossing χ=A.This indicates that δCMis no longer a singlevalued function with respect to χ.Therefore,correction is needed.

    Figure 22 has the same plasma discharge parameters and toroidal magnetic field settings as figure 21.The only difference lies in the initial polarization parameters (χ0,ψ0) set asThe code demonstrates that,when the CM effect is sufficiently strong,it has the ability to reverse the handedness of a wave that is initially launched with circular polarization.A noticeable alignment between the corrected results and the definition is evident,leading to a reduction in error from 48.34 % (whenx0=-0.225 m) to approximately 0.5 %.Whenx0>0,the sign consistency of ψ before and after passing through the plasma ensures that it does not enter the domain of equation (64),thus no correction is needed.

    Figure 22.Correction of computational equation.The results calculated using equation (5) are represented by dashed lines.Due to the consistency of the plasma discharge parameters and the toroidal magnetic field between the scenarios in figure 21,the ‘Definition’curves in the two plots are identical.The result of the calculation solely using equation (55) is represented as ‘ ×’ and the result of the correction using both equation (55) and equation (64) is represented as ‘+’.

    Although the nonlinear expression is precise,its limitation lies in the complex square root operation required for signal demodulation.Moreover,ensuring the consistency of the sign of ψ during the measurement process is challenging,and neglecting the Faraday rotation (FR) effect will lead to a significant increase in the calculation error of χ.This makes it difficult to use as a diagnostic tool for real-time feedback of electron density during operation and control of fusion reactors.In the balance between accuracy and real-time capability,practical plasma discharges often do not adopt ellipsometry measurements.

    5.Conclusion

    This study combines theoretical derivation and simulations to systematically investigate the feasibility of measuring electron density using the Cotton-Mouton (CM) effect on EAST.For laser wavelength of 432.6 μm,fringe jumps due to phase shift over 2 π can be avoided to provide reliable density measurement for machine operation.The study clarifies that elliptical modulation enables the direct measurement of the phase variation of the second components2of the Stokes parameters,providing information aboutThe measurement principle for heterodyne detection in the case of small CM effects is presented,along with the design of an optical setup.In a laboratory environment,the power of the two orthogonal laser beams received by the detector may not be strictly equal.It is shown that non-equal power radiation only affects the DC offset of the light intensity,and does not affect the phase variation of measured signals.It is further confirmed that the higher-order term errors introduced by the Faraday rotation effect on CM effect are small (about 1 %)for current parameters of EAST.However,in the case of large CM effects (δCM~O(1)),the higher-order terms of the Mueller matrix will introduce significant errors.

    Current plasma density diagnostics based on polarimeter interferometers consider the ellipticity parameters as linear functions of phase shifts.As future fusion devices improve in their parameters,this approach will introduce increasingly large errors.This paper reveals the coupling between the ellipticity parameters χ and ψ,the potential errors they may introduce in plasma density diagnostics based on the CM effect.The study provides an explicit expression for the CM effect phase shift δCMas a function of χ,decoupled from the azimuthal angle ψ,when the Faraday rotation effect can be neglected.

    In summary,for EAST experiment parameters,a verticalview Cotton-Mouton polarimeter interferometer is capable of providing reliable density measurement with excellent time resolution (a few MHz).

    Acknowledgments

    This work was financially supported by National Natural Science Foundation of China (No.12127809).

    Appendix A

    The following section will demonstrate the analytical expression for the first-order expansion ofs2(z) whenW~O(1).Whens(z=0) is modulated with a significant CM effect,numerical integration becomes inconvenient.

    In such cases,the zeroth-order Mueller matrix differential equation is introduced to facilitate the analysis [20]

    Given t he initial conditionM(0)(z=0)=E,whereEdenotes the identity matrix,the integration yields:

    Based on the reference [17],it can be derived that the first-order terms of Mueller matrix do not affect the phase ofs2(z).

    Similarly to equation (26),the second components2of the Stokes parameters is given by

    Namely,for even values ofn,M(n) has the same distribution of zeros asM(0)and,for odd values ofn,the same asM(1).It can be observed that the non-zero elementsM22andM23of the second-order matrixM(2)will affect the phase ofs2(z)[14].

    Appendix B

    Now,the indefinite integral of equation (55) is to be calculated.From equation (54),

    The positive sign is taken when sin2χ>0 and the negative sign is taken when sin2χ<0.

    Derivation of equation (77)

    substituting equation (79) into equation (76),

    derivation of equation (81) yields

    substituting equation (84) into equation (80),

    Also because

    When χ0,ψ0are known,the constantCcan be found.The sign of δCMdepends on the positivity or negativity of sin2χ,which is taken to be positive when sin2χ>0 and positive when sin2χ<0 is taken as negative.

    猜你喜歡
    劉海衛(wèi)星
    以小見大 以情動(dòng)人
    miniSAR遙感衛(wèi)星
    如何確定衛(wèi)星的位置?
    軍事文摘(2021年16期)2021-11-05 08:48:58
    Establish a Three-dimensional Fluorescent Fingerprint Database of Traditional Chinese Medicines
    Hubbard model on an anisotropic checkerboard lattice at finite temperatures:Magnetic and metal–insulator transitions
    學(xué)生天地(2019年36期)2019-08-25 08:59:52
    靜止衛(wèi)星派
    科學(xué)家(2019年3期)2019-08-18 09:47:43
    Puma" suede shoes with a focus on the Product variables
    只靠劉海就能實(shí)現(xiàn)的超簡(jiǎn)單變身方法!
    競(jìng)射導(dǎo)航衛(wèi)星為哪般
    太空探索(2015年6期)2015-07-12 12:48:29
    欧美精品一区二区免费开放| 国产黄色视频一区二区在线观看| 国产91av在线免费观看| 久久这里有精品视频免费| 色综合色国产| 亚洲精品久久久久久婷婷小说| 婷婷色综合大香蕉| 成人漫画全彩无遮挡| 六月丁香七月| 韩国高清视频一区二区三区| 联通29元200g的流量卡| 久久久久久久久久久免费av| 亚洲精品亚洲一区二区| 国产精品成人在线| av福利片在线观看| 亚洲欧美日韩东京热| 国产精品一区二区三区四区免费观看| av在线老鸭窝| 国产成人aa在线观看| 欧美zozozo另类| 欧美精品国产亚洲| 看非洲黑人一级黄片| .国产精品久久| 国产精品.久久久| 亚洲av中文字字幕乱码综合| 男人狂女人下面高潮的视频| 久久99热6这里只有精品| 天美传媒精品一区二区| 亚洲怡红院男人天堂| 男女免费视频国产| 免费观看av网站的网址| 一本—道久久a久久精品蜜桃钙片| 久久精品久久精品一区二区三区| 国产老妇伦熟女老妇高清| 亚洲美女搞黄在线观看| 国产在线免费精品| 五月伊人婷婷丁香| 久久久久久久精品精品| 国产一级毛片在线| 一个人看视频在线观看www免费| 亚洲av国产av综合av卡| 亚洲欧美日韩无卡精品| 身体一侧抽搐| kizo精华| 久久久久久久大尺度免费视频| 国产精品一区二区在线观看99| 久久久久精品性色| 一边亲一边摸免费视频| 汤姆久久久久久久影院中文字幕| 五月天丁香电影| 亚洲av日韩在线播放| 夫妻性生交免费视频一级片| 久久久久国产网址| 中文在线观看免费www的网站| 精品一区二区免费观看| 狠狠精品人妻久久久久久综合| 久久6这里有精品| 国内少妇人妻偷人精品xxx网站| 国产欧美日韩精品一区二区| 男人和女人高潮做爰伦理| 日韩制服骚丝袜av| 永久免费av网站大全| 九九在线视频观看精品| 91狼人影院| 亚洲图色成人| 一级二级三级毛片免费看| 亚洲成人中文字幕在线播放| 三级国产精品欧美在线观看| 欧美变态另类bdsm刘玥| 中文字幕人妻熟人妻熟丝袜美| 免费观看的影片在线观看| 久久 成人 亚洲| 国产男女内射视频| 人妻一区二区av| 亚洲精品中文字幕在线视频 | 亚洲在久久综合| 日本免费在线观看一区| 久久精品熟女亚洲av麻豆精品| 亚洲欧美日韩另类电影网站 | 黄色视频在线播放观看不卡| 国产成人午夜福利电影在线观看| av国产免费在线观看| 国产高清国产精品国产三级 | 免费高清在线观看视频在线观看| 在线观看av片永久免费下载| 精品人妻视频免费看| 啦啦啦在线观看免费高清www| av国产精品久久久久影院| 18禁裸乳无遮挡免费网站照片| 极品教师在线视频| 国产午夜精品一二区理论片| 99视频精品全部免费 在线| 久久99蜜桃精品久久| 国产永久视频网站| 日本黄色日本黄色录像| 美女cb高潮喷水在线观看| 男的添女的下面高潮视频| 色网站视频免费| 黄色配什么色好看| 久久亚洲国产成人精品v| 午夜视频国产福利| 蜜桃久久精品国产亚洲av| 日本黄色片子视频| 男人和女人高潮做爰伦理| 狂野欧美激情性bbbbbb| av女优亚洲男人天堂| 久热久热在线精品观看| 欧美+日韩+精品| 久久人妻熟女aⅴ| 高清黄色对白视频在线免费看 | 国产精品国产三级国产av玫瑰| 日韩三级伦理在线观看| 亚洲国产精品专区欧美| 草草在线视频免费看| av在线蜜桃| 国产高清不卡午夜福利| 国产在线免费精品| 高清午夜精品一区二区三区| 国产成人91sexporn| 日本猛色少妇xxxxx猛交久久| 国模一区二区三区四区视频| 丝瓜视频免费看黄片| 熟女av电影| 熟女电影av网| 国产精品成人在线| 国产精品精品国产色婷婷| av不卡在线播放| 国产v大片淫在线免费观看| 国产日韩欧美亚洲二区| 性色av一级| 51国产日韩欧美| 国产高清三级在线| 人体艺术视频欧美日本| 亚洲av男天堂| 欧美成人一区二区免费高清观看| 久久精品国产亚洲av天美| 日韩欧美一区视频在线观看 | 国产精品熟女久久久久浪| 特大巨黑吊av在线直播| 尾随美女入室| 日韩av在线免费看完整版不卡| 亚洲精品乱久久久久久| av国产久精品久网站免费入址| 麻豆精品久久久久久蜜桃| 成年女人在线观看亚洲视频| 欧美xxⅹ黑人| 免费播放大片免费观看视频在线观看| 久久久久久久精品精品| 国产无遮挡羞羞视频在线观看| 午夜免费观看性视频| 久久久久久久久久久免费av| 久久99热这里只有精品18| 久久久久久久精品精品| 特大巨黑吊av在线直播| 国产色爽女视频免费观看| 日韩欧美一区视频在线观看 | 婷婷色综合www| 18禁在线无遮挡免费观看视频| 久久影院123| 我要看日韩黄色一级片| 日日撸夜夜添| 国产综合精华液| 成年人午夜在线观看视频| 亚洲成人一二三区av| 人人妻人人看人人澡| 天堂中文最新版在线下载| 国产精品一及| 在线观看国产h片| 欧美 日韩 精品 国产| 尾随美女入室| 国产极品天堂在线| 亚洲精品国产色婷婷电影| 久久久久精品性色| 少妇的逼水好多| 久久久亚洲精品成人影院| 亚洲国产欧美在线一区| 性色avwww在线观看| 中文字幕亚洲精品专区| 亚洲欧美日韩东京热| 欧美精品一区二区免费开放| 一级毛片久久久久久久久女| 国产v大片淫在线免费观看| 又黄又爽又刺激的免费视频.| 老司机影院毛片| 国产精品偷伦视频观看了| 免费久久久久久久精品成人欧美视频 | 一区二区三区乱码不卡18| 中文字幕亚洲精品专区| 成人亚洲精品一区在线观看 | 免费不卡的大黄色大毛片视频在线观看| 嫩草影院新地址| 青春草国产在线视频| 久热这里只有精品99| 日本黄大片高清| 黑人高潮一二区| 欧美性感艳星| 免费少妇av软件| 免费观看av网站的网址| 免费观看性生交大片5| 涩涩av久久男人的天堂| 亚洲欧美一区二区三区黑人 | 国产亚洲午夜精品一区二区久久| 秋霞在线观看毛片| 成人高潮视频无遮挡免费网站| 99久久中文字幕三级久久日本| 精品久久久精品久久久| av黄色大香蕉| 男的添女的下面高潮视频| 一区二区av电影网| 一本色道久久久久久精品综合| 亚洲国产欧美在线一区| 搡女人真爽免费视频火全软件| 一级毛片 在线播放| 亚洲av男天堂| 各种免费的搞黄视频| 各种免费的搞黄视频| 高清不卡的av网站| 亚洲图色成人| 一级毛片黄色毛片免费观看视频| 各种免费的搞黄视频| 下体分泌物呈黄色| 一级毛片电影观看| 中文乱码字字幕精品一区二区三区| 国产成人aa在线观看| av不卡在线播放| 国产片特级美女逼逼视频| a级毛色黄片| 联通29元200g的流量卡| 国产黄频视频在线观看| 日韩免费高清中文字幕av| 国产成人精品一,二区| 女人久久www免费人成看片| 久久久久视频综合| 少妇人妻精品综合一区二区| 国产高清不卡午夜福利| 精品久久久久久久久av| 国产中年淑女户外野战色| 国产白丝娇喘喷水9色精品| 国内揄拍国产精品人妻在线| 日本一二三区视频观看| 国产一区二区三区综合在线观看 | 日本黄色片子视频| 大香蕉久久网| 亚洲欧美成人综合另类久久久| 六月丁香七月| 久久精品熟女亚洲av麻豆精品| 久久精品国产自在天天线| 久久韩国三级中文字幕| 久久精品国产亚洲网站| 九草在线视频观看| 亚洲国产精品专区欧美| 亚洲欧美一区二区三区国产| 日韩电影二区| 最近最新中文字幕免费大全7| 日韩亚洲欧美综合| 看十八女毛片水多多多| 18禁裸乳无遮挡动漫免费视频| 国产片特级美女逼逼视频| 亚洲av二区三区四区| 爱豆传媒免费全集在线观看| 人妻系列 视频| 欧美区成人在线视频| 国产精品国产三级国产av玫瑰| 香蕉精品网在线| 亚洲色图av天堂| 国产黄色免费在线视频| 中文字幕精品免费在线观看视频 | 久久久久性生活片| 伦理电影大哥的女人| 日韩制服骚丝袜av| 91在线精品国自产拍蜜月| 一级av片app| 老司机影院毛片| 啦啦啦在线观看免费高清www| 国产中年淑女户外野战色| 精品少妇黑人巨大在线播放| 91aial.com中文字幕在线观看| 久久久久网色| 伦理电影大哥的女人| 久久ye,这里只有精品| 少妇 在线观看| 大话2 男鬼变身卡| 亚洲国产最新在线播放| 久久女婷五月综合色啪小说| 午夜福利视频精品| 亚洲无线观看免费| 少妇的逼好多水| 亚洲av免费高清在线观看| 亚洲av成人精品一区久久| 91久久精品电影网| 久久久精品94久久精品| 91精品国产九色| 日韩电影二区| 午夜福利在线在线| 久久ye,这里只有精品| 日本黄色片子视频| 亚洲成人av在线免费| 国产欧美另类精品又又久久亚洲欧美| 国产精品嫩草影院av在线观看| 国产在线一区二区三区精| 黑丝袜美女国产一区| 我的女老师完整版在线观看| 成人无遮挡网站| 新久久久久国产一级毛片| 亚洲成人av在线免费| 三级国产精品欧美在线观看| 亚洲av中文av极速乱| 欧美日韩精品成人综合77777| 国产精品免费大片| 久久ye,这里只有精品| 一级毛片我不卡| 精品一区二区三卡| 免费人妻精品一区二区三区视频| 91精品国产国语对白视频| 直男gayav资源| 人体艺术视频欧美日本| 亚洲精品国产av蜜桃| 久久精品国产鲁丝片午夜精品| 国产久久久一区二区三区| 麻豆成人av视频| 一个人看的www免费观看视频| 菩萨蛮人人尽说江南好唐韦庄| 国语对白做爰xxxⅹ性视频网站| 人妻 亚洲 视频| 五月伊人婷婷丁香| 中文乱码字字幕精品一区二区三区| 国产成人freesex在线| 亚洲av.av天堂| 免费观看无遮挡的男女| 一级av片app| av在线观看视频网站免费| 国产免费视频播放在线视频| 日韩成人伦理影院| 久久99热这里只有精品18| a 毛片基地| 亚洲国产日韩一区二区| 少妇人妻久久综合中文| 97精品久久久久久久久久精品| 亚洲欧美日韩东京热| 一级毛片 在线播放| 插阴视频在线观看视频| 欧美成人a在线观看| 国产在线一区二区三区精| 国产高清三级在线| 一本久久精品| 日韩中文字幕视频在线看片 | 国内揄拍国产精品人妻在线| 91久久精品国产一区二区三区| 一级毛片我不卡| 亚洲精品中文字幕在线视频 | 晚上一个人看的免费电影| 18+在线观看网站| 男女啪啪激烈高潮av片| 男女国产视频网站| 一级黄片播放器| 国产伦精品一区二区三区四那| 亚洲,欧美,日韩| 欧美区成人在线视频| 深夜a级毛片| 国产男女内射视频| 成年人午夜在线观看视频| 少妇 在线观看| 国产日韩欧美亚洲二区| 国产欧美亚洲国产| 国产午夜精品久久久久久一区二区三区| 亚洲av中文av极速乱| 亚洲综合色惰| 国产淫片久久久久久久久| 亚洲国产精品国产精品| 美女中出高潮动态图| 国产黄色视频一区二区在线观看| 成人综合一区亚洲| 国产精品一区二区在线观看99| 亚洲精品中文字幕在线视频 | av播播在线观看一区| 十分钟在线观看高清视频www | 九草在线视频观看| 国产在线视频一区二区| 大话2 男鬼变身卡| 91久久精品电影网| 永久网站在线| 久久久久精品久久久久真实原创| 国产精品女同一区二区软件| 简卡轻食公司| 国产 一区 欧美 日韩| 人体艺术视频欧美日本| 日韩av在线免费看完整版不卡| 久久人人爽人人爽人人片va| 男女啪啪激烈高潮av片| 毛片女人毛片| 亚洲欧美清纯卡通| 少妇人妻精品综合一区二区| 色婷婷久久久亚洲欧美| 晚上一个人看的免费电影| 欧美精品亚洲一区二区| 亚洲欧美成人综合另类久久久| 夜夜看夜夜爽夜夜摸| 免费av不卡在线播放| 午夜福利在线在线| 亚洲国产av新网站| 激情五月婷婷亚洲| 蜜桃亚洲精品一区二区三区| 高清视频免费观看一区二区| www.色视频.com| 国产一级毛片在线| 亚洲欧美日韩东京热| 九色成人免费人妻av| 久久精品久久久久久久性| 亚洲国产日韩一区二区| 日韩国内少妇激情av| 亚洲天堂av无毛| 青春草亚洲视频在线观看| 夜夜骑夜夜射夜夜干| 国产精品嫩草影院av在线观看| 国产成人freesex在线| 久久久成人免费电影| 秋霞在线观看毛片| 亚洲国产日韩一区二区| 夫妻午夜视频| 国产探花极品一区二区| 黑人猛操日本美女一级片| 伊人久久国产一区二区| 五月伊人婷婷丁香| 免费在线观看成人毛片| 亚洲av欧美aⅴ国产| 狂野欧美激情性xxxx在线观看| 国国产精品蜜臀av免费| 久久婷婷青草| 综合色丁香网| 美女高潮的动态| 国产成人一区二区在线| 老熟女久久久| 国产精品国产av在线观看| 亚洲精品视频女| 一级毛片电影观看| 日韩一区二区三区影片| 少妇的逼水好多| 看十八女毛片水多多多| 国产真实伦视频高清在线观看| 亚洲aⅴ乱码一区二区在线播放| 亚洲精品国产av蜜桃| 高清午夜精品一区二区三区| 国产男人的电影天堂91| 18禁裸乳无遮挡免费网站照片| 国产亚洲一区二区精品| 欧美日本视频| 欧美激情国产日韩精品一区| 亚洲av日韩在线播放| 黑人高潮一二区| 亚洲精品国产色婷婷电影| 热99国产精品久久久久久7| 在线观看人妻少妇| 亚洲欧美一区二区三区黑人 | 色哟哟·www| 啦啦啦视频在线资源免费观看| 欧美日韩在线观看h| 91精品伊人久久大香线蕉| 国产在线男女| 又粗又硬又长又爽又黄的视频| 香蕉精品网在线| 女人久久www免费人成看片| 国产伦理片在线播放av一区| xxx大片免费视频| 成人美女网站在线观看视频| 亚洲国产高清在线一区二区三| 大香蕉久久网| 国产黄频视频在线观看| 在线观看美女被高潮喷水网站| 香蕉精品网在线| 在线观看免费视频网站a站| 亚洲av国产av综合av卡| 搡老乐熟女国产| 97精品久久久久久久久久精品| 看免费成人av毛片| 亚洲怡红院男人天堂| 大香蕉97超碰在线| 91精品国产九色| 免费av不卡在线播放| 国产精品偷伦视频观看了| 91久久精品国产一区二区成人| 国产老妇伦熟女老妇高清| 极品少妇高潮喷水抽搐| 直男gayav资源| 26uuu在线亚洲综合色| 亚洲人与动物交配视频| 五月开心婷婷网| 九色成人免费人妻av| 亚洲欧洲日产国产| 男女边吃奶边做爰视频| 中文字幕亚洲精品专区| 男的添女的下面高潮视频| 又粗又硬又长又爽又黄的视频| 久久精品久久久久久久性| 男人添女人高潮全过程视频| www.色视频.com| 国产黄片视频在线免费观看| 国产精品麻豆人妻色哟哟久久| 亚洲欧美精品专区久久| 亚洲成色77777| 国产精品一区二区在线不卡| 最近的中文字幕免费完整| 蜜桃在线观看..| 我的老师免费观看完整版| 久久人人爽人人爽人人片va| 韩国av在线不卡| 精品国产露脸久久av麻豆| 亚洲av电影在线观看一区二区三区| 中文字幕免费在线视频6| 99视频精品全部免费 在线| 美女国产视频在线观看| 日本一二三区视频观看| 免费少妇av软件| 汤姆久久久久久久影院中文字幕| 久久综合国产亚洲精品| 欧美丝袜亚洲另类| 亚洲久久久国产精品| 国产精品爽爽va在线观看网站| 黄色欧美视频在线观看| 亚洲精品久久久久久婷婷小说| 国产亚洲一区二区精品| 国产男女超爽视频在线观看| videossex国产| 一个人看的www免费观看视频| 一区二区三区精品91| 最近中文字幕高清免费大全6| 欧美人与善性xxx| 国产精品一区二区在线不卡| 人妻 亚洲 视频| 国产精品三级大全| 日韩欧美一区视频在线观看 | 午夜精品国产一区二区电影| av免费观看日本| 青春草国产在线视频| 婷婷色麻豆天堂久久| 亚洲国产精品成人久久小说| 一个人看的www免费观看视频| 亚洲欧美一区二区三区国产| 91久久精品电影网| 亚洲av免费高清在线观看| 欧美最新免费一区二区三区| 色吧在线观看| 日本欧美视频一区| 天美传媒精品一区二区| 身体一侧抽搐| 欧美少妇被猛烈插入视频| 日韩一区二区视频免费看| 夫妻午夜视频| 久久久久久久国产电影| 日本午夜av视频| 美女cb高潮喷水在线观看| 国产熟女欧美一区二区| 99久久精品热视频| 丝袜喷水一区| 日本-黄色视频高清免费观看| 夜夜骑夜夜射夜夜干| 人人妻人人爽人人添夜夜欢视频 | 久久久久久伊人网av| 内地一区二区视频在线| 国产精品无大码| 男人狂女人下面高潮的视频| 国产精品成人在线| 热99国产精品久久久久久7| 美女主播在线视频| 中文字幕免费在线视频6| 五月伊人婷婷丁香| 精品人妻一区二区三区麻豆| 日韩大片免费观看网站| 欧美日韩国产mv在线观看视频 | 99久久中文字幕三级久久日本| 久久久久久伊人网av| 91精品国产九色| 日本与韩国留学比较| 我要看黄色一级片免费的| 十八禁网站网址无遮挡 | 亚洲一级一片aⅴ在线观看| 热re99久久精品国产66热6| 亚洲伊人久久精品综合| 国产大屁股一区二区在线视频| 美女脱内裤让男人舔精品视频| 国语对白做爰xxxⅹ性视频网站| 久久鲁丝午夜福利片| 偷拍熟女少妇极品色| 亚洲av中文字字幕乱码综合| 亚洲av福利一区| 国产欧美日韩一区二区三区在线 | 少妇人妻一区二区三区视频| 亚洲国产精品国产精品| 一本—道久久a久久精品蜜桃钙片| 特大巨黑吊av在线直播| 日本vs欧美在线观看视频 | 久久精品国产亚洲网站| 日本色播在线视频| 97精品久久久久久久久久精品| 国产真实伦视频高清在线观看| 亚洲精品国产av蜜桃| 国产 一区精品| 亚洲精品一区蜜桃| 婷婷色综合大香蕉| 高清在线视频一区二区三区| 五月伊人婷婷丁香| 内地一区二区视频在线| 日日摸夜夜添夜夜添av毛片| 黑人高潮一二区| 久久精品国产鲁丝片午夜精品| 亚洲无线观看免费| 成人高潮视频无遮挡免费网站| 欧美日韩一区二区视频在线观看视频在线| 99热这里只有精品一区| 国产在线视频一区二区| 亚洲天堂av无毛| 在线免费观看不下载黄p国产| 久久国产精品大桥未久av | 亚洲精品久久午夜乱码| 麻豆精品久久久久久蜜桃| 身体一侧抽搐|